A dynamical system method for solving the split convex feasibility problem

In this paper a dynamical system model is proposed for solving the split convex feasibility problem. Under mild conditions, it is shown that the proposed dynamical system globally converges to a solution of the split convex feasibility problem. An exponential convergence is obtained provided that the bounded linear regularity property is satisfied. The validity and transient behavior of the dynamical system is demonstrated by several numerical examples. The method proposed in this paper can be regarded as not only a continuous version but also an interior version of the known \begin{document}$ CQ $\end{document} -method for solving the split convex feasibility problem.

1. Introduction. The split convex feasibility problem (SCF P ) due to Censor and Elfving [16] is to find x ∈ R n such that x ∈ C and Ax ∈ Q, where C ⊂ R n and Q ⊂ R m are two nonempty closed convex sets, and A ⊂ R m×n is a real matrix. SCF P arises in many applications, such as the phase retrieval problems [16], signal processing [13], image reconstruction [28,35], the intensitymodulated radiation therapy [15], and the gene regulatory network inference arising in systems biology [44]. A more general model is the following prototypical split inverse problem (SIP ) introduced by Censor et al. [17]: where X and Y are two spaces, A : X → Y is a bounded linear operator, IP 1 is an inverse problem formulated in the space X and IP 2 is another inverse problem formulated in the space Y . SCF P is the first instance of a SIP in which two inverse problems are convex feasible problems. In fact, many models of inverse problems can be cast in SIP by choosing different inverse problems for IP 1 and IP 2 . For details, we refer the reader to [17]. In the literature, many numerical algorithms have been proposed for solving SCF P (see e.g., [41,38,18,2,44,19,37,25]). Among them, a classical algorithm for SCF P is the known CQ-algorithm proposed by Byrne [13,14] which is defined by: for given x 0 , where γ is a suitable positive constant, A T is the transpose of A, I is the identical operator, P C and P Q are the orthogonal projections onto C and Q, respectively. Byrne's CQ-algorithm is an effective and popular method for solving SCF P because it involves only easily calculated orthogonal projections and it requires no matrix inversions. It is worth mentioning that the CQ-algorithm is actually a projected gradient method applied to the following constrained minimization problem: which encapsulates the split convex feasibility problem (SCF P ). Among the existing works on CQ-algorithm and its variations for SCF P , a great deal of attention has been paid to the convergence of algorithms, but only a few papers deal with the convergence rate. Under the bounded linear regularity condition, Wang et al. [44] studied the linear convergence rate of the varying stepsize CQ-algorithm for SCF P . Qu et al. [37] proposed a Newton projection method (which can also be reviewed as a variation of CQ-algorithm) for SCF P and analyzed the convergence properties of the method. Especially, under some assumptions which are necessary for the case where the projection operator is nonsmooth, Qu et al. [37] proved that the sequence of iterates generated by the algorithm superlinearly converges to a regular solution of SCF P . On the other hand, dynamical system approaches were proposed by scholars for solving optimization problems, Wardropian user equilibria, variational inequality problems and inclusion problems. Pyne [36] first proposed the dynamical system approach for solving the linear programming problem on an electronic analogue computer. Xia and Wang [47,46] presented neural networks for solving linear and quadratic programming problems. Attouch et al. [3] proposed a gradientprojection dynamical system for solving constrained optimization problems. Effati et al. [21] proposed a projection neural network for solving convex programming problems. Friesz et al. [24] presented a dynamical system model for calculating static Wardropian user equilibria on congested networks with fully general demand and cost structure. Xia and Wang [48] proposed a projection dynamical system for solving the variational inequality problem. Attouch and Svaiter [5] proposed a continuous dynamical Newton-like approach for solving the monotone inclusion problem. Abbas et al. [1] proposed a dynamical system for solving the inclusion problem with the sum of a convex subdifferential and a monotone cocoercive operator. For more results on dynamical system approaches, we refer the reader to [4,10,23,26,32,34,33,49,51] and the references therein. Motivated by the above works, in this paper, we attempt to propose a dynamical system model for solving SCF P . In order to present our model, we rewrite CQalgorithm (2) as By bringing stepsize h n into the above equation, we get Set (3), we propose the following dynamical system for solving SCF P : where λ > 0 is a scaling factor. By the above arguments, system (4) can be regarded as a continuous version of CQ-algorithm (2). Comparison with the CQ-algorithm, system (4) enjoys advantages and much stronger properties. In general, the sequence generated by the CQ-algorithm (2) is not in the interior int C of C even if the initial point x 0 ∈ int C. We will see that the trajectory of system (4) with the initial point being an interior always lies in the interior of C; see Remark 4. So system (4) can be considered as not only a continuous version of the CQ-algorithm (2) but also an interior version. As known, the explicit discretizing scheme for system (4) may lead to the CQ-algorithm. In fact, we can also obtain other numerical algorithms for (SCF P ) by using different discretizing schemes. For example, taking a fixed time step size h > 0 and setting t n = nh, x n = x(t n ), the central difference scheme for system (4) gives: Furthermore, system (4) has the advantage that it can be solved by different methods such as the finite difference methods and the finite element method. This will been seen in Section 5.
The remainder of this paper is organized as follows. In Section 2, we recall some concepts and preliminary results. In Section 3, we investigate SCF P by proving that system (4) globally converges to a solution of SCF P . We further obtain an exponential convergence under the bounded linear regularity condition. In Section 4, we apply results of Section 3 to the linear equation problem. In Section 5, numerical examples are reported to illustrate the effectiveness and performance of the dynamical system approach. Finally, we end the paper with a conclusion.

2.
Preliminaries. In this section, we recall some concepts and results which will be used in the sequel.
and T is called nonexpansive if it is 1-Lipschitz continuous. (b) T is called r-inverse strongly monotone on R n if there exists r > 0 such that T is called firmly nonexpansive if it is 1-inverse strongly monotone.
(c) T : R n −→ R n is called an α-averaged mapping if it can be written as where α ∈ (0, 1), I denotes the identity mapping, and S : R n −→ R n is a nonexpansive mapping.

Remark 1. (i) T is firmly nonexpansive if and only if
(ii) T is nonexpansive if it is averaged or firmly nonexpansive. (iii) T is firmly nonexpansive if and only if I − T is firmly nonexpansive. (iv) P C and I − P C are both firmly nonexpansive, i.e., and , then U : R n → R n is an averaged mapping.
Lemma 2.4. (See [14,38]) Suppose that the solution set Λ of SCF P is nonempty. Let x * ∈ R n . Then the following statements are equivalent: (i) x * solves the SCF P .
Definition 2.5. (See [44]) Suppose the solution set Λ of SCF P is nonempty. We say that SCF P satisfies the bounded linear regularity property if, for any r > 0 such that Λ ∩ B(0, r) = ∅, there exists κ r > 0 such that where B(x, r) and B(x, r) denote the open and closed ball with center at x and radius r, respectively.
Remark that regularity-type conditions have been widely used to analyze the convergence rates of many algorithms; see [7,12,30,44] and references therein. Definition 2.6. (See [42,46]) Let Φ : R n → R n . Consider the following dynamical system (7) is called stable if, for any ε > 0, there exists δ > 0 such that for every x 0 ∈ B(x * , δ), the solution x(t) of system (7) with initial point x(t 0 ) = x 0 exists and is contained in denotes the open ball with center at x * and radius ε.
(c) A stable equilibrium point x * of (7) is called asymptotically stable if there exists δ > 0 such that, for every solution (7) is called globally convergent to a set Ω ⊂ R n if, for any solution x(t) of (7), Remark 2. By (a) of Definition 2.6 and Lemma 2.4, the solution set Λ of SCF P coincides with the set of all the equilibrium points for system (4). [40,42]) An equilibrium point x * of system (7) is stable if and only if there exists a Lyapunov function V (x) : R n → R for the dynamical system (7) at x * , which is defined as: for some neighborhood B(x * ) of x * , V (x) satisfies: has continuous first order partial derivatives; (iii) its time derivative along any state trajectory x(t) of dynamical system (7) is is locally Lipschitz continuous at x 0 then the solution is unique, and if Φ(x) is Lipschitz continuous on R n then τ can be extended to ∞.
Lemma 2.9. (See [27]) Let t 0 ∈ R and h : (t 0 , +∞) → R n be a uniformly continuous function, e.g., for every real > 0 there exists > 0 such that for every To end this section, we recall the following identity will be used in the paper (see Corollary 2.14 of [8]): 3. Convergence and stability. In this section, we shall analyze convergence of a solution of dynamical system (7), and then derive some convergence rate. Before that, we first prove the existence and uniqueness of solution. In what follows, we always suppose that γ ∈ (0, 2/ A 2 ).
Proof. Set T = P C • U − I. Then (4) can be written as By Lemma 2.2 and (ii) of Remark 1, U is nonexpansive. We claim that T is Lipschitz continuous. Indeed, for any x, y ∈ R n , if follows that By Lemma 2.8, there exists a unique solution x(t) of (4) satisfying In the following theorem, we discuss the stability of system (4) in the sense of Lyapunov and the convergence of (I − P Q )Ax(t) and Theorem 3.2. For arbitrary initial point x 0 ∈ R n , let x(t) be the unique solution of dynamical system (4) with x(t 0 ) = x 0 . Then, for every equilibrium point x * ∈ Λ, the following statements hold: Since x * ∈ Λ, by Lemma 2.4 we have For any initial point x 0 ∈ R n , by Theorem 3.1, there exists a unique solution x(t) of system (4) with x(t 0 ) = x 0 . We claim that dV (x(t)) dt ≤ 0. Indeed, from (4), (10) and (11) we have It follows from (5), (6), (10) and (12) that ≤ 0.
By Lemma 2.7, the equilibrium point x * of system (4) is stable in the sense of Lyapunov. Next, we verify the statements (ii) and (iii). Since V (x(t)) is nonincreasing as t → +∞ and bounded (it is nonnegative), we deduce that lim t→+∞ V (x(t)) exists and the trajectory x(t) is bounded. Integrate (13) from t 0 to +∞ to get So (ii) and (iii) hold. Now we proof (iv). Integrating (9) from t 1 to t 2 and using mean value theorem, there existst ∈ [t 1 , t 2 ] such that This along with the bounded of x(t) and the Lipschitz property of T yield the assertion (iv). And finally we shall show that (v) and (vi) are true. It is clear that I − P Q is Lipschitz continuous and A is a bounded linear operator. Combining with (iv) we have that (I − P Q )Ax is Lipschitz continuous. Therefore, the assertions (v) and (vi) are true by Lemma 2.9.
In Theorem 3.2, we have studied the stability of system (4). In fact, we can further obtain the global convergence of system (4) and a convergence rate Theorem 3.3. Let x(t) be the solution of system (4) with x(t 0 ) = x 0 ∈ R n . Then the following statements are true: (i) x(t) converges to some equilibrium point of system (4) as t → ∞.
. To prove (i), we divide the proof into two steps.
Step 1: We shall show that x−x * is nonincreasing on [0, +∞), for every x * ∈ Λ. Set By the classical derivation chain rule, we have where the fifth equality is from (8). So x − x * is nonincreasing on [0, +∞), for every x * ∈ Λ.
Step 2. We shall show that lim t→+∞ dist(x(t), Λ) = 0. It is worth noting that x(t) − x * is nonincreasing and x(t) is bounded on [0, +∞). Letx be a cluster of x(t) and t n → ∞ such that x(t n ) →x. From (vi) in Theorem 3.2, we have P C • Ux =x. By Lemma 2.4,x ∈ Λ. By replacing x * withx in (14), we get that lim t→+∞ x(t) −x exists. So x(t) →x as t → +∞.
In order to prove (ii), we first shall show that P C •Ux(t)−x(t) is nonincreasing on [t 0 , +∞). By the Rademacher's theorem (see Theorem 3.16 in [22]) and the differentiability of x(t), we know that P C • Ux(t) is differentiable almost everywhere with respect to t. Since P C • U is nonexpansive, by (b) of Remark 1 of [11], we get that holds almost everywhere. Set G(x(t)) = 1 2 P C • U(x(t)) − x(t) 2 . Deriving G(x(t)) with respect to time, from (4) and (16) we have that holds almost everywhere. This together with the continuity of P C • Ux(t) implies that P C • Ux(t) − x(t) is nonincreasing on [0, +∞).
Then using the monotonicity of G(x(t)) and the conclusion (vi) in Theorem 3.2, we obtain → 0 ( as t → ∞).
As a consequence, Remark 3. An abstract dynamical system governed by a non-expansive mapping T instead of P C •U has been already considered in Bot and Csetnek [11]. As pointed out by one reviewer that: (a) (iii) and (vi) of Theorem 3.2, and Theorem 3.3 can follow directly from Theorem 6 and Theorem 10 of Bot and Csetnek [11](although our method is different from the one in [11]); (b) The dynamics proposed in [10] for a structured convex minimization problem is also able to solve the split convex feasibility problem.
It is worth mentioning that Bot and Csetnek [11] did not consider the exponential convergence of the dynamical system. Next, we further derive the exponential convergence of the trajectory of system (4) with an initial point x 0 ∈ C when the bounded linear regularity property holds. To do this, we need the following lemma. Lemma 3.4. Let f : [0, ∞) → R n be a continuous mapping such that f (t) ∈ C for all t ≥ t 0 and let x(t) be a trajectory of the following dynamical system: where λ > 0 and C is a nonempty, closed and convex set. Then we have x(t) ∈ C for all t ≥ t 0 .
Proof. From (17) we get Integrating it from t 0 to t, we have This yields We claim that Indeed, let [t i , t i+1 ], i = 0, 1, · · · , N be any tagged partition of [t 0 , t], := max{ i : i = 1, 2, · · · , N } be the mesh of such a tagged partition with i being the width of sub-interval [t i , t i+1 ]. Let τ i ∈ [t i , t i+1 ]. By the definition of integral we get This together with (18) implies that x(t) ∈ C for all t ≥ t 0 .
Remark 4. The well known methods for solving approximatively system (4) are finite difference methods. Taking time step size t n > 0, the explicit difference scheme for system (4) gives: For λ t n ≡ 1, the CQ-algorithm (2) is recovered in (19). Comparison with the CQalgorithm (2), system (4) enjoys much stronger property. In general, the sequence generated by the CQ-algorithm (2) is not in the interior int C (respectively, relative interior ri C) of C even if the initial point x 0 ∈ int C (respectively, ri C) while system (4) is an interior method as soon as the initial point is strictly feasible. Indeed, by the argument (18) of Lemma 3.4, we have x(t) ∈ int C (respectively, ri C) for all t ≥ t 0 if x 0 ∈ int C (respectively, ri C). In fact, the discrete approximation (19) with λ t n < 1 also enjoys this interior property. Now, with Lemma 3.4 in hand, we derive the exponential convergence of the trajectory of system (4).
Theorem 3.5. Let x(t) be the solution of system (4) with x(t 0 ) = x 0 ∈ C and x ∈ Λ such that x(t) →x as t → ∞(This is guaranteed by Theorem 3.3). Suppose that SCF P satisfies the bounded linear regularity property. Then there existt ≥ t 0 and β 0 > 0 such that Proof. Set f (t) = P C • Ux(t). Clearly, f is continuous and f (t) ∈ C for all t ∈ t 0 . By Lemma 3.4, x(t) ∈ C for all t ≥ t 0 . Define By Lemma 2.4, we have P C • U(P Λ x(t)) = P Λ x(t). Deriving W (x(t)) with respect to time t, we obtain where the first equality follows from (a) of Lemma 2.3. Recalling that U = I − γA T (I − P Q )A, we get We now estimate the last two terms of (21). By (iii) of Lemma 2.4 and the definition of L, we get that From (6) and (iv) of Lemma 2.4 we have that Substituting these two inequalities into (21), we have This together with (20) yields that Since x(t) ∈ C for all t ≥ t 0 , x(t) →x as t → ∞ and SCF P satisfies the bounded linear regularity property, there exist κ > 0 andt ≥ t 0 such that It follows from (22) and (23) that It follows that which implies thatW (t) is nonincreasing on [t, ∞). From (25) and (26) we get This yields By (22), dist 2 (x(t), Λ) is nonincreasing on [t 0 , ∞). This together with the last inequality implies that Remark 5. In some sense, Theorem 3.5 can be viewed as a continuous form of Theorem 2.3 of [44] where the linear convergence of a varying stepsize CQ algorithm was established under the bounded linear regularity property.

4.
Application to the linear equation problem. When C = R n and Q = {b}, the split feasible problem (1) reduces to the following linear equation problem (LEP): find x ∈ R n such that In this case, Ax), and CQ-algorithm algorithm (2) reduces to the following Landweber method [31]: for x 0 arbitrary. As a consequence, system (4) reduces to the following dynamical system: Naturally, the Landweber method (28) can be regarded as an explicit discrete-time form of system (29). When A is a matrix with full column rank, an implicit discretetime scheme of system (29) yields the Piccard algorithm for (27): where ρ > 0. As an application of results in Section 3, we solve LEP (27) via system (29). In what follows, we always suppose that the solution set E of (27) is nonempty.
Proof. It is a direct consequence of Theorem 3.1.
Theorem 4.2. For arbitrary initial point x 0 ∈ R n , let x(t) be the unique solution of dynamical system (29) with x(t 0 ) = x0. Then the following the following statements hold: (i) for every equilibrium point x * ∈ E of system (29) is stable; Proof. It is a direct consequence of Theorem 3.2. Theorem 4.3. Let x(t) be the solution of system (29) with x(t 0 ) = x 0 ∈ R n . Then the following statements are true: (i) x(t) converges to some equilibrium point of system (29) as t → ∞.
Proof. Part (i) is a direct consequence of (i) of Theorem 3.3. Part (ii) follows from (ii) of Theorem 3.3 and the following fact: Theorem 4.4. Let x(t) be the solution of system (29) with x(t 0 ) = x 0 ∈ R n and x ∈ E such that x(t) →x as t → ∞(This is guaranteed by Theorem 4.3). Then there exist t ≥ t 0 and β 0 > 0 such that ∀t ≥t.
Example 2. (See [45, Example 1]) We consider SCF P (1) as following: the matrix A = (a i,j ) 2×2 and a i,j ∈ (0, 1) are generated randomly, the nonempty closed convex Take γ = 1 A 2 and λ = 1 3 . Figure   Here, the explicit difference method is exploited to solve approximately dynamical system (4). By means of algebraic calculation, it is easy to verify that x * = [9, −1, −6] T is a unique solution for Example 3. Set γ = 1 A 2 and λ = 35 that satisfy the operating conditions. Figure 6, where the explicit difference method is exploited, displays that the trajectories of dynamical system (4) with 20 random initial points all globally converge to x * . Figure 7 and Figure 8 depict the asymptotical behavior the trajectory of dynamical system (4) with initial points x 0 = [−3, 5, 2] T where the finite element method with λ = 1 and the Piccard algorithm with ρ = 1 5 A 2 are exploited to solve dynamical system (4), respectively.
Example 4. In this example we use our dynamical system approach for solving the problem of recovering a noisy sparse signal from a limited number of sampling which was considered in [25,Subsection 5.1]. Let x ∈ R n be a K−sparse signal and K << n. The sampling matrix A ∈ R m×n (m < n) is stimulated by standard Gaussian distribution and vector b = Ax + e with e being an additional noisy. The task is to recover the signal x from the data b. This problem can be modeled as the following LASSO problem [43]: where t > 0 is a given constant. Taking C = {x : x 1 ≤ t} and Q = {b}, the above LASSO problem is transformed into SCF P (1). So we can solve the problem via dynamical system (4).
In the implementation of dynamical system (4) we choose the following set of parameters: the initial point x 0 ∈ R n of dynamical system (4) and the matrix A ∈ R m×n (m < n) are generated randomly. The true signal x ∈ R n with K non−zero elements is generated from uniform distribution in the interval [−2, 2] and vector b = Ax. Take m = 2 12 , n = 2 13 , t = K = 50, γ = 1.9 A 2 and λ = 5. Figure 9 describes the asymptotical behavior of dynamical system (4) where ode45 is exploited to solve (4). It is shown in Figure 9 that the trajectory x(t) of dynamic system (4) converges. Figure 10 depicts the recovery of the true signal by the proposed dynamical system approach. Figure 11 displays the objective function value against the time for the LASSO problem solved through the dynamical system (4) with different choices of λ.
6. Conclusion. In this paper we propose a projection dynamical system for solving SCF P . We demonstrate that the trajectory of the projection dynamical system is globally convergent to a solution of SCF P . It is also shown that the fixed point residual enjoys the o( 1 √ t ) convergence rate. Under the bounded linear regularity condition, we further obtain an exponential convergence for the trajectory of the the projection dynamical system. The method proposed in this paper can be regarded as not only a continuous version but also an interior version of the known CQmethod for solving the split convex feasibility problem.