- Email: [email protected]
Projected randomized Kaczmarz methods
Journal Pre-proof Projected randomized Kaczmarz methods Nianci Wu, Hua Xiang
PII: DOI: Reference:
S0377-0427(19)30677-6 https://doi.org/10.1016/j.cam.2019.112672 CAM 112672
To appear in:
Journal of Computational and Applied Mathematics
Received date : 12 March 2019 Revised date : 10 August 2019 Please cite this article as: N. Wu and H. Xiang, Projected randomized Kaczmarz methods, Journal of Computational and Applied Mathematics (2020), doi: https://doi.org/10.1016/j.cam.2019.112672. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Elsevier B.V. All rights reserved.
Journal Pre-proof
PROJECTED RANDOMIZED KACZMARZ METHODS NIANCI WU
†
∗
AND HUA XIANG †
of
Abstract. For a constrained linear ill-posed problems, the projected algebra reconstruction technique (PART) usually has a faster convergence rate than the projected simultaneous iterative reconstruction technique (PSIRT), but can not be applied to the inconsistent case. In this work we propose the projected randomized Kaczmarz method (PRK) to further accelerate the convergence of PART and analyze its convergence. With the noise-free right-hand side, we show that the PRK method can converge exponentially in expectation to the solution for the consistent system. With the noisy right-hand side, we present some insights into the semiconvergence property of the PRK method and propose an extended version of PRK (PREK) method to seek the least squares solution. We prove that the PREK method can converge exponentially in expectation to the least squares solution for the inconsistent systems. Key words. projected algebra reconstruction technique, projected randomized iterative method, exponential convergence, semiconvergence
p ro
AMS subject classifications. 15A23, 65F30, 65F22, 65F10, 65R32
1. Introduction. Consider the constrained discrete ill-posed problems Ax = b subject to x ∈ C,
(1.1)
(1.2)
urn
al
Pr e-
where A ∈ Rm×n (m ≥ n), b ∈ Rm , and C is the nonempty closed convex set. The discrete ill-posed problems are characterized by the coefficient matrices with a very large condition number. This implies that the naive solution is very sensitive to any perturbation of the right-hand side, representing the errors in the data. The numerical solution of these systems requires the use of iterative regularization methods. Furthermore, iterative methods are usually more suitable than direct methods for a large scale problem. Some numerical techniques in the real-life applications can be found in different fields (see [6, 9, 30, 36, 40, 42]). It is often necessary to incorporate additional constrains to the iterative regularization methods, which can typically lead to smaller iteration errors in computational tomography as well as many other image problems (see [1, 14, 16, 17]). The common closed convex set includes nonnegativity constrain, box constrain, finite support constrain and bandlimiting constrain (see [37]). Assuming that S and T are two closed subspaces of a Hilbert space H, to find the best approximation to x from S ∩ T , the method of alternating projection (MAP) [19] suggests that we can first project x onto S, the obtained element is then projected onto T , and the algorithm projects alternatively onto S and T . In this way it generates a sequence of elements that converges to PS∩T (x). Combining the projection methods with the reconstruction techniques, e.g., simultaneous iterative reconstruction technique (SIRT) [14], in each step of the MAP method, Elfving, Hansen and Nikazad [16] give the following projected SIRT (PSIRT) method, which has extensive applications in image restoration, matrix rank minimization and matrix completion (see [16, 29, 39, 43, 46]). Algorithm PSIRT. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 21 = xk−1 + λk−1 AT M (b − Axk−1 ),
xk = PC (xk− 21 ),
Jo
where {λk }∞ k=0 is a sequence of positive relaxation parameters, M is a symmetric positive definite matrix and PC is the metric projection onto the nonempty closed convex set C. Some well-known methods can be written in the form (1.2) for appropriate choices of the matrix M . The following three choices of M are quite popular: (1) M = Im where Im is the identity Pn matrix of order m, 1 (2) M = m diag(1/kai k2 ) where aTi denotes the i-th row of A, and (3) M = diag(1/ j=1 Nj a2ij ) where Nj is the number of nonzeros in the j-th column of A and aij refers to the (i, j) entry of A. These three cases correspond to the classical projected Landweber method [27, 29, 43, 46], the projected Cimmino method [10, 16], and the projected component averaging (CAV) method [1, 8, 16], respectively. Let h be the spectral radius of the matrix AT M A and assume that 0 ≤ ² ≤ λk ≤ (2 − ²)/h. If ² > 0, then the iterates of SIRT with and without projections converge to a solution of min kAx − bkM [14], where kuk2M = uT M u is the norm induced by the symmetric positive definite matrix M for any vector u. Moreover, developing appropriate step size strategies can significantly improve the convergence rate of SIRT, for more details, refer to [14, 32]. The Kaczmarz method [26] proposed in 1937 is a simple yet powerful technique for solving the linear system. In tomography research it is also known as the algebra reconstruction technique (ART) [24]. ART ∗ This work was supported by the National Natural Science Foundation of China under Grant No. 11571265 and NSFC-RGC No. 11661161017. † School of Mathematics and Statistics, Wuhan University, Wuhan 430072, PR China. ([email protected]).
1
Journal Pre-proof 2
Projected randomized Kaczmarz method
is a particular case of row projection methods that suits for parallel computations and large scale problems since each step only requires one row of matrix A and no matrix-vector products [7], while SIRT makes 1 use of all the equations simultaneously. For example, by taking M = m diag(1/kai k2 ) and the relaxation parameter as constant 1, SIRT automatically reduces to the Cimmino method, which computes the next iterative by averaging all the projections of the previous iterates. The iterate schemes of the Cimmino and ART methods are given by
of
xk− 21
¶ m m µ 1 X 1 X bi − aTi xk−1 = PHi (xk−1 ) ≡ xk−1 + ai m i=1 m i=1 kai k2 ¶ m µ 1 X bi − aTi xk−1 a = xk−1 + AT M (b − Axk−1 ) = xk−1 + i m i=1 kai k2
and bi − aTi xk−1 ai kai k2
(i = 1, 2, · · · , m, 1, 2, · · · , m, 1, 2, · · · ),
p ro
xk− 12 = PHi (xk−1 ) ≡ xk−1 +
respectively, where PHi is the metric projection onto the hyperplane Hi = {x ∈ Rn : hai , xi = bi },
Pr e-
and bi is the i-th component of the right-hand side vector b ∈ Rm . In addition, we find that if we have B = AAT at the beginning, then ART only costs 4n + 1 flopping operations (flops), while the Cimmino method requires for 4mn + m flops, at each iteration step. In the recent work of Elfving, Hansen and Nikazad [17], they have found empirically that ART usually has a faster convergence rate than some SIRT such as the Landweber and Cimmino methods. The geometry of the Cimmino and ART methods are also shown in Figure 1.1 when m = 2, which can intuitively explain the fast convergence of the ART method compared to the Cimmino method. H1
xk-1
xk-1
al
H1
urn
H2
(a) The Cimmino method
H2
(b) The ART method
Fig. 1.1. The geometry of the Cimmino and ART methods. The hyperplane Hi = {x ∈ Rn : hai , xi = bi } for i = 1, 2.
(1.3)
Jo
By substituting SIRT with ART in the PSIRT method, Andersen and Hansen present the following projected ART (PART) method in [1]. Algorithm PART. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 12 = PHi (xk−1 ),
xk = PC (xk− 12 ),
i = 1, 2, · · · , m, 1, 2, · · · , m, 1, 2, · · · ,
where PHi and PC are the metric projections onto the hyperplane Hi and the nonempty closed convex set C, respectively. The theoretically convergence rate of ART is very challenging, which is reflected by the fact that the convergence rate of the method depends strongly on the ordering of the equations. Strohmer and Vershynin in an influential paper [38] use the rows of A in ART with a random order, called the randomized Kaczmarz (RK) method below, rather than the given order. This strategy can guarantee an exponential convergence rate1 of the RK method for consistent linear systems, where the convergence rate depends on the condition number. This result is then extended and refined in various directions [2, 3, 4, 5, 21, 28, 31, 41, 45]. 1 A randomized iterative method that computes a sequence of random vectors {x }∞ k k=0 is said to converge in expectation to a vector x∗ if and only if E[kxk − x∗ k2 ] → 0 as k → ∞, where the expectation is taken over the random choice of the methods.
Journal Pre-proof 3
NIANCI WU AND HUA XIANG
of
Inspired by the PART and RK methods, we propose the projected randomized Kaczmarz (PRK) method to solve the constrained linear ill-posed problems. First, the PRK method is given in section 2, including its convergence analyses for the noise-free and noisy right-hand sides. Second, we present an extended version of the PRK (PREK) method that can converge exponentially in expectation to the least square solution for the constrained inconsistent systems in section 3. Finally, we demonstrate the convergence behaviors of the proposed algorithms with the numerical examples from medical tomography. Throughout the paper k·k denotes the vector and matrix 2-norm, and k·kF denotes the matrix Frobenius norm. For a rectangular matrix A, we denote by A† the Moore-Penrose pseudoinverse of A [20], R(A) the range of A and N (AT ) the null space of AT . Given any vector b ∈ Rm , we can uniquely write it as b = br + bn , where br = AA† b and bn = (Im − AA† )b are the projections of b onto R(A) and N (AT ), respectively. We indicate by Ek−1 [·] the expected conditional on the k-th iteration, and from the law of the iterated expectation we have E[Ek−1 [·]] = E[·]. Contributions. We comment on the main contributions of this work as follows.
Pr e-
p ro
(i) A new scheme. The PART is a popular way to solve the constrained linear system, which is under the framework of the MAP method. This method suggests that at the first-half step, we use ART to the linear system Ax = b, and at the second-half step, the obtained element is projected onto the set C. Then this process continues alternatively. Based on the fact that the RK method has a faster convergence rate than ART, we replace ART with RK in the PART method and get the PRK method to solve the constrained linear system. We show that the PRK method possesses the expected exponential convergence rate. (ii) Semiconvergence analysis. We apply the PRK method to the case where the constrained linear system is inconsistent, which has the noisy right-hand side (i.e., bn 6= 0). We observe that the PRK method will not converge to the least squares (LS) solution of the inconsistent system. Some insights into its semiconvergence property are presented also. To investigate how the noisy part bn is propagated during the PRK iterations, we split the total error into iteration error and noise error, and give an upper estimate of the noise error. (iii) A extended variant for the LS solution. We utilize the randomized orthogonal projection method to find the approximation of bn and employ the PRK method to the approximate consistent system Ax ≈ b−bn . Then the PREK method is obtained, which is a new LS solver of the inconsistent constrained linear system. We prove that the PREK method can converge exponentially in expectation to the LS solution for the inconsistent system. 2. The projected randomized Kaczmarz method. In this section, we first present the PRK method to solve the constrained consistent linear ill-posed problems. Then, its convergence analyses are followed for the noise-free and noisy right-hand sides separately.
(2.1)
urn
al
2.1. The PRK method. Motivated by the PART process, the PRK method is mainly based on the next two generic step procedures. We first project the iterate vector xk−1 onto the i-th hyperplane Hi in a random order, and then project the current vector onto the closed convex set C. Continuing with this process alternatively, the PRK method can be formulated as follows. Algorithm PRK. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 21 = PHi (xk−1 ),
xk = PC (xk− 12 ),
Jo
where PHi and PC are the metric projections onto the hyperplane Hi and the nonempty closed convex set C, respectively, the row index i is chosen from the set [m] := {1, 2, · · · , m} at random with probability pi = kai k2 /kAk2F . Remark 1. Let Sk = {x ∈ Rn : Aτk x = bτk , τk ⊂ [m]}, where bτk ∈ Rτk , Aτk ∈ Rτk ×n is formed by the rows of A, then for any y ∈ / Sk the projection onto Sk is given by PSk (y) = y + ATτk (Aτk ATτk )† (bτk − Aτk y).
Instead of just using one row aTi of A in (2.1), we can work on several rows simultaneously, it is natural to extend the PRK method to block variants. Different from ART, which selects the row index cyclically in a deterministic fashion, the RK method selects the row index at random with a probability criterion. Thus, the PRK method can be viewed as the random version of the PART method. Additionally, at the first-half iteration step, the PART and PRK methods need the same flops. From the geometry viewpoint, we know that if we project the current vector xk− 12 onto the hyperplane Hi (i = 1, 2, · · · , m) as far as possible at the first-half step of the PART method, then it could lead to a faster convergence rate. This also can be seen from the Figure 1.1(b). In designing a random version of the PART method, i.e., the PRK method, we aim at grasping larger entries of the residual vector and selecting the working row according to a probability criterion, rather than in their given order. It is possible to avoid a hyperplane, which has a small distance from xk−1 , to be selected.
Journal Pre-proof 4
Projected randomized Kaczmarz method
k I 2 E[keIk k2 ] ≤ (1 − κ−2 A ) · ke0 k ,
p ro
(2.2)
of
2.2. Convergence analysis of PRK with noise-free right-hand side. Let x∗ be the solution of (1.1) and xk denote the k-th iterate of Algorithm PRK using the noise-free right-hand side, i.e., b = Ax∗ , then bi = aTi x∗ for i = 1, 2, · · · , m. Define the iteration error eIk := xk − x∗ , which is also known as data error, see [14, 16, 17]. Next, we will give an upper estimate on E[keIk k2 ] for the PRK method. This represents our first main theoretical result. Lemma 2.1. [19, Proposition 4.1] If C is a closed and convex set in any Hilbert space H, then the projection operator PC : H → C is non-expansive, i.e., for any elements α, β ∈ H, kPC (α) − PC (β)k ≤ kα − βk. By utilizing the above lemma, we can deduce the convergence property of the PRK method with noisefree right-hand side and establish the following theorem. Theorem 2.2. Let the constrained linear system (1.1), with the coefficient matrix A ∈ Rm×n and the right-hand side b ∈ Rm , be consistent. Then from any initial guess x0 ∈ R(AT ), the iteration sequence ∗ {xk }∞ k=0 generated by the PRK method converges to the solution x in expectation. Moreover, the mean squared iteration error for k = 0, 1, 2, · · · satisfies
where the generalized condition number κA = kAkF · kA† k. Proof. Based on the non-expansivity of the projection operator in Lemma 2.1, it holds that m m X X kai k2 kai k2 ∗ 2 ∗ 2 = · kPC [PHi (xk−1 )] − PC (x )k ≤ 2 2 · kPHi (xk−1 ) − x k kAk kAk F F i=1 i=1
= (2.3) =
m X kai k2 bi − aTi xk−1 ∗ · k(x − x ) + ai k2 k−1 2 2 kAk ka k i F i=1
Pr e-
Ek−1 [keIk k2 ]
m X kai k2 ai aTi · k(In − )eIk−1 k2 2 2 kAk ka k i F i=1
m X kai k2 ai aTi AAT I I I T · (eIk−1 )T (In − )e = )e . = (e ) (I − n k−1 k−1 2 kAkF kai k2 kAk2F k−1 i=1
al
By induction, xk−1 is in the row space of A for any k when x0 ∈ R(AT ). In addition, the same is true for x∗ by the definition of psedoinverse. This means that eIk−1 = xk−1 − x∗ is also in the row space of A. Therefore, we can have µ ¶ σr2 (A) AAT I I 2 I T )e ≤ 1− · keIk−1 k2 = (1 − κ−2 (ek−1 ) (In − A ) · kek−1 k . kAk2F k−1 kAk2F
urn
By taking the full expectation for both sides of the inequality (2.3), we obtain I 2 E[keIk k2 ] ≤ (1 − κ−2 A ) · E[kek−1 k ].
Jo
Repeating the same argument k − 1 times we get the estimate (2.2). Strohmer and Vershynin in [38] prove a similar exponential upper bound on the expected convergence rate for the noise-free right-hand side about the RK method. Theorem 2.2 indicates that the convergence rate of the PRK method is not less than that of the RK method due to the non-expansivity of the projection operator. Besides, the expected convergence rate of the RK method in [38] is 1−b κ−2 , where κ b = kAkF ·kA−1 k is the scaled condition number introduced by Demmel [11]. It is based on the restrictive assumption that A has full column rank, while the convergence rate of the PRK method only depends on the generalized condition number κA without imposing this assumption. In practice, one usually observes the typical semiconvergence phenomenon for iterative regularization methods, e.g., unprojected/projected SIRT in [16, 14], ART in [1, 17] and the general iteration scheme in [15]. That is, the iterative procedure of these methods are damaged by noise. The iteration vector approaches a regularized solution in the early stage, while it is deteriorated by noise in the later iterations. Next we will study the semiconvergence behavior of the PRK method. 2.3. Semiconvergence analysis of PRK with noisy right-hand side. If the constrained linear system (1.1) is inconsistent, we use x∗ = A† b ∈ C to represent its LS solution. Then it holds that the ek denote the k-th iterates of the PRK method using the noisy right-hand side b = Ax∗ + bn . Let xk and x ek − x∗ and the noise error noise-free and noisy right-hand sides, respectively. Define the total error e ek := x N ek − xk , where the noise error is also known as approximation error, see [14, 16, 17]. ek := x
Journal Pre-proof 5
NIANCI WU AND HUA XIANG
I Clearly, the total error e ek satisfies that e ek = eN k + ek . Classical convergence theorem usually focuses on eIk (e.g., Theorem 2.2), while assuming that eN is negligible. However, this is not the case for ill-posed and k noisy problems. The interplay of these two error terms determines the final total error. We next investigate how the noisy part bn in the right-hand side is propagated during the PRK iterations. Theorem 2.3. Let the linear system (1.1) be with the coefficient matrix A ∈ Rm×n . The iteration sequences {xk }∞ xk }∞ k=0 and {e k=0 are generated by the PRK method with noise-free right-hand side b = m br ∈ R and noisy right-hand side b = br + bn ∈ Rm , respectively. Assuming that the initial guesses e0 = x0 ∈ R(AT ), then the mean squared noise error eN ek − xk is bounded by x k =x
√
k kbn k2 · , σr (A) kAkF
where σr (A) is the smallest singular value of A Proof. Direct calculation yields that µ
ek−1 + x
ek−1 eTi (br + bn ) − aTi x ai 2 kai k
¶
µ
¶ eTi bn − aTi xk−1 a . i kai k2
p ro
eN k = PC
of
2 E[keN k k ]≤
(2.4)
− PC
xk−1 +
Using again the non-expansivity of the projection operator, we can have
m m X X kai k2 kai k2 ai aTi eTi bn N 2 N · ke · k(I − k ≤ )e + a k2 n k k−1 2 2 2 2 i kAk kAk ka k ka k i i F F i=1 i=1 · ¸ m 2 X ai aTi (eTi bn )2 kai k N T N · (e ) (I − )e + = n k−1 kAk2F kai k2 k−1 kai k2 i=1 "m # m X kai k2 X ai aTi (eTi bn )2 T N = (eN · (I − ) e + n k−1 ) k−1 kAk2F kai k2 kAk2F i=1 i=1
Pr e-
2 Ek−1 [keN k k ]=
N 2 ≤ (1 − κ−2 A ) · kek−1 k +
kbn k2 , kAk2F
where ei is the i-th column of identity matrix Im with order m. By taking the full expectation for both sides of the above inequality, it indicates that k−1 2 X kbn k2 −2 k −2 j kbn k N 2 ≤ (1 − κ ) · ke k + . (1 − κ ) · 0 A A kAk2F kAk2F j=0
al
−2 2 N 2 E[keN k k ] ≤ (1 − κA ) · E[kek−1 k ] +
urn
Assuming x = 1/κA ∈ (0, 1), the above inequality reads 2 E[keN k k ]≤
k−1 X j=0
j (1 − κ−2 A ) ·
e0 = x0 . Let y = 1 − x2 ∈ (0, 1), then since x
Jo
(2.5)
1 − (1 − x2 )k = x
s
(1 − y k )2 ≤ 1−y
kbn k2 1 − (1 − x2 )k kbn k2 = · κ · , A kAk2F x kAk2F
s
p √ 1 − yk = 1 + y + · · · + y k−1 ≤ k, 1−y
and the bound (2.4) comes out. This completes the proof. √ Note that the bound in Theorem 2.3 includes a factor k, which is similar to the formula for the propagation of the noise error in ART and SIRT, see [14, 17], and is valid also for their projected versions, see [16, 17]. The iteration number k and the smallest singular value σr (A) of A play the role of the regularized parameters. The bound (2.4) might be quite pessimistic since it holds for any closed convex set C and the estimates used in (2.5) are rather crude. As stated in Theorem 2.3, when the system (1.1) is perturbed by noise or no longer consistent, it is a pity that the PRK method will not converge to the LS solution of the inconsistent system. It is because that at the k-th iterate the PRK method projects the current vector xk−1 onto the hyperplane aTi x = eTi b = eTi (br +bn ) rather than aTi x = eTi br for i = 1, · · · , m. To resolve this problem, we provide the following variant of the PRK method.
Journal Pre-proof 6
Projected randomized Kaczmarz method
zk = zk−1 −
(3.1)
Aj ATj zk−1 , kAj k2
of
3. The projected randomized extended Kaczmarz method. In this section we will present a projected randomized iterative LS solver that converges in expectation to the LS solution of the inconsistent linear system with constraints. For the general inconsistent linear system, Zouzias and Freris in [45] modify the RK method and obtain a random version of extended Kaczmarz method motivated by the work of Popa in [33, 34]. We call it REKZF. One step of the REK-ZF method consists of two stages. The first stage uses one step RK update for the linear system Ax = b − zk−1 and generates xk . The second stage applies one step randomized orthogonal projection update to the homogeneous linear system AT z = 0 and outputs zk , which admits the following iterate scheme
Pr e-
p ro
where Aj is the j-th column of A and the index j is chosen at random. Actually, the iterate vector zk is the projection of zk−1 onto the hyperplane ATj z = 0, and it is equivalent to apply one RK update to AT z = 0. The REK-ZF method has an expected exponential convergence rate and immediately results in many works devoted to various aspects of the new algorithm, which includes distinct bounds on their performance and its extensions. Du presents a slightly different REK-ZF (denoted by REK-Du) method in [12], which generates xk by one-step RK update for the linear system Ax = b − zk (used in REK-Du) instead of Ax = b − zk−1 (used in REK-ZF) from xk−1 , because zk is a better approximation of bn than zk−1 . Recently, Bai and Wu [5] compute zk from (3.1) in a given cyclic order rather than in a dynamic random order, which can obtain the approximation of bn more accurately, and present the partial REK-ZF method. We call it REK-BW. It is proved that the REK-Du and REK-BW methods all have the smaller upper bounds for the expected convergence rate than that of the REK-ZF method (see [12, Thm.2] and [5, Thm.3.1]). 3.1. The PREK method. We know from Theorem 2.3 that the reason why the PRK method fails to converge to LS solution of the inconsistent problems is that bn 6= 0. Then, we split the PREK method into two computational stages. The first stage outputs zk by utilizing one step randomized orthogonal projection update to the linear system AT z = 0. The second stage uses one step PRK update to the approximate consistent linear system Ax = b − zk with x ∈ C. Note that if the nonempty closed convex set C = Rn , the PREK method could reduce to the REK-Du method. The pseudo-code of the PREK method is given in the following algorithm. Algorithm PREK. Given two initial guesses x0 ∈ Rm and z0 ∈ b + R(A), for k = 1, 2, · · · until ∞ {xk }k=0 converges, compute zk = PHe j (zk−1 ),
xk− 12 = PHb i (xk−1 ),
al
(3.2)
xk = PC (xk− 21 ),
where PHe j , PHb i and PC are the metric projections onto the hyperplanes
urn
e j = {z ∈ Rm : hAj , zi = 0}, H
b i = {x ∈ Rn : hai , xi = eT (b − zk )} H i
and the nonempty closed convex set C, respectively, the column index j is chosen from the set [n] at random with probability qj = kAj k2 /kAk2F and the row index i is chosen from the set [m] at random with probability pi = kai k2 /kAk2F .
(3.3)
Jo
3.2. Convergence analysis of PREK. For the convergence property of the PREK method, we can establish the following theorem. First, we review a known result needed in the remaining parts of this section. Lemma 3.1. [12, Thm.1] Let A ∈ Rm×n and b ∈ Rm . Let zk denote the k-th iterate of the randomized orthogonal projection method applied to AT z = 0 with initial guess z0 ∈ b + R(A). The iteration sequence {zk }∞ k=0 converges to the solution br in expectation. Moreover, the mean squared error for k = 0, 1, 2, · · · satisfies E[kzk − bn k2 ] ≤ ρk kz0 − bn k2 ,
where ρ = 1 − σr2 (A)/kAk2F = 1 − κ−2 A . Theorem 3.2. Let the constrained linear system (1.1), with the coefficient matrix A ∈ Rm×n and the right-hand side b = br + bn ∈ Rm , be inconsistent. Then from any initial guesses x0 ∈ R(AT ) and z0 ∈ b + R(A), the iteration sequence {xk }∞ k=0 generated by the PREK method converges to the solution x∗ = A† b in expectation. Moreover, the mean squared error for k = 0, 1, 2, · · · satisfies (3.4)
E[kxk − x∗ k2 ] ≤
kρk kz0 − bn k2 + ρk kx0 − x∗ k2 . kAk2F
Journal Pre-proof 7
NIANCI WU AND HUA XIANG
Proof. After the k-th iteration of the PREK method, it yields that kxk − x∗ k2 = kPC (xk−1/2 ) − PC (x∗ )k2
eTi (b − zk ) − aTi xk−1 ai − A† bk2 kai k2 eT (bn − zk ) eT br − aTi xk−1 = kxk−1 + i ai − A† b + i ai k2 2 kai k kai k2 eT AA† b − aTi xk−1 eT (bn − zk ) = kxk−1 + i ai − A† b + i ai k2 2 kai k kai k2 ¶ µ eT (bn − zk ) ai aTi (xk−1 − A† b) + i = k In − a i k2 2 kai k kai k2 ¶ µ eTi (bn − zk ) ai aTi † 2 = k In − (x − A b)k + k a i k2 , k−1 kai k2 kai k2 ≤ kxk−1/2 − x∗ k2 = kxk−1 +
p ro
of
(3.5)
(3.6)
Pr e-
where the second line comes from the the non-expansivity of the projection operator, the´third line follows ³ ai aT T from the fact that b = bn + br and the last line is from the orthogonality ai In − kai ki2 = 0. Like [12], we use Ek−1 [·] to denote the expected value conditional on the first k−1 iterations of the PREK method, i.e., Ek−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 ], where j` and i` (` = 1, · · · , k−1) mean that the `-th column and row are chosen as Ejk−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 , ik ] and Eik−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 , jk ], respectively. Then, by the law of total expectation, i.e., Ek−1 [·] = Ejk−1 [Eik−1 [·]], it follows that ¸ ¸¸ · · ¸¸ · T · T · T e (bn − zk ) (ei (bn − zk ))2 ei (bn − zk ) j j 2 2 i i a k = E a k = E E E k Ek−1 k i i i k−1 k−1 k−1 k−1 kai k2 kai k2 kai k2 # "m · ¸ X kai k2 (eT (bn − zk ))2 kbn − zk k2 j j i · = Ek−1 = Ek−1 kAk2F kai k2 kAk2F i=1 =
£ ¤ 1 ρk 2 E kb − z k ≤ kbn − z0 k2 , k−1 n k kAk2F kAk2F
al
where the last line is from Lemma 3.1. Similar to the derivation of Theorem 2.2, we have that (3.7) · µ ¶ ¸ £ ¤ £ ¤ ai aTi † 2 † 2 (xk−1 − A b)k ≤ (1 − κ−2 = ρ Ek−1 kxk−1 − A† bk2 . Ek−1 k In − A ) · Ek−1 kxk−1 − A bk 2 kai k
urn
In accordance with formulas (3.6) and (3.7), the formula (3.5) then further results in the estimate · T ¸ · µ ¶ ¸ £ ¤ e (bn − zk ) ai aTi 2 † 2 Ek−1 kxk − x∗ k2 ≤ Ek−1 k i a k + E k I − (x − A b)k i k−1 n k−1 kai k2 kai k2 £ ¤ ρk ≤ kbn − z0 k2 + ρ Ek−1 kxk−1 − A† bk2 . 2 kAkF By taking the full expectation on both sides of the above inequality, it follows that
Jo
£ ¤ E kxk − x∗ k2 ≤
£ ¤ ρk kbn − z0 k2 + ρ E kxk−1 − A† bk2 . kAk2F
Then we immediately obtain the estimate in (3.4) by induction on the iteration index k. 4. Numerical experiments. In this section we investigate some numerical tests taken from the field of tomography image reconstruction. To discretize the tomography problems, we should divide the square domain into a grid of N × N pixels numbered from 1 to n (n = N 2 ). The j-th cell is assigned a constant value xj of the attenuation or travel time for j = 1, 2, · · · , n, such that the unknown vector x is a discrete version of the sought function, i.e., x = [x1 , · · · , xn ]T . Define the coefficient matrix A = [aij ] ∈ Rm×n , where aij is the length of the i-th ray through the j-th cell. Let the number of rays be m, then the measurement bi of the attenuation or travel time for the i-th ray can be written as bi =
n X j=1
aij xj = aTi x for i = 1, 2, · · · , m.
Journal Pre-proof 8
Projected randomized Kaczmarz method
For more details, refer to [22]. As the set C we use the positive orthant such that we enforce nonnegativity on the constructions. The orthogonal projector PC onto the positive orthant can be defined as follows. ½ yi , yi ≥ 0, [PC (y)]i = 0, yi < 0,
Pr e-
p ro
of
for any vector y = [y1 , · · · , yn ]T . Thus, a constrained discrete ill-posed problem is determined. Consider the Sheep-Logan phantom, its exact solution is shown in Figure 4.1(a). Based on various scanning modes, e.g., the fan-beam, spherical radon and parallel-beam method, where their geometries are shown with a few of rays in Figure 4.1((b)-(d)), the MATLAB package (AIR Tools, available from https://github.com/jakobsj/AIRToolsII) in [22, 23] presents three 2-dimension (2D) tomography test examples. They are fancurvedtomo, sphericaltomo and paralleltomo. These examples generate a sparse matrix A, an exact solution x∗ and an exact right-hand side b = Ax∗ . It gives us many flexibilities to adjust the input parameters, such as the size of the discrete domain (N ), the number of rays (p), the projection angles in degrees (θ) and the number of sources (s), where m = ps and n = N 2 . Therefore, we can assign values to these parameters and obtain various specific instantiations about the constrained linear systems. In this work, we typically set N = 40, p = 40, θ = 0 : 2 : 358 and s = 180. The resulting matrix is of size 7200 × 1600. All computations are started from the initial guess x0 = 0, and terminated once the relative solution error at the k-th iterate satisfies kxk − x∗ k ≤ kx∗ k · 10−6 or the number of iteration steps exceeds 105 . The numerical tests are performed on a Founder desktop PC with Intel(R) Core(TM) i5-7500 CPU 3.40 GHz by MATLAB R2016(a) with a machine precision of 10−16 .
(a) Sheep-Logan Phantom
(b) fan-beam
(c) spherical radon
(d) parallel-beam
al
Fig. 4.1. The exact solution of the Sheep-Logan phantom (a). The geometries of fan-beam (b), spherical radon (c) and parallel-beam (d), with a few of rays.
1
1
10
10
ART RK PART PRK Error
Error
100
ART RK PART PRK
100
urn
Error
ART RK PART PRK
-1
10
10-2
100
10-1
-3
10
2
4 6 Iteration number
Jo
(a) fancurvedtomo
8
10
4
×10
2
4 6 Iteration number
(b) sphericaltomo
8
10 4
×10
1
2
3
4 5 6 7 Iteration number
8
9
10 4
×10
(c) paralleltomo
Fig. 4.2. Error histories of the ART, RK, PART and PRK methods for fancurvedtomo(a), sphericaltomo (b) and paralleltomo (c) with noise-free right-hand side.
Case 1: noise-free right-hand side. In this case, we utilize the ART, RK, PART and PRK methods to solve the constrained consistent linear system. The errors kxk − x∗ k versus iteration number k for both methods with noise-free right-hand side are shown in Figure 4.2. In this figure, we can observe that the PRK method converges faster than the ART, PART and RK methods. Case 2: noisy right-hand side. To examine the semiconvergence phenomenon of the PRK method, like [17, 25], the noisy right-hand side b is defined as b = A · x∗ + η · kbk · ξ/kξk, where η is the relative noise level and ξ is a random vector generated by the MATLAB function randn. In this case the system (1.1) is no longer consistent, we first depict the convergence behavior of the ART, RK, PART and PRK methods in Figure 4.3 – Figure 4.5 with three different relative noise levels η = 1 × 10−2 , 5 × 10−2 and 1 × 10−1 , which
Journal Pre-proof 9
NIANCI WU AND HUA XIANG
9 8 7
ART RK PART PRK
9
ART RK PART PRK
6
7 6 Error
Error
5 Error
ART RK PART PRK
8
4
5
3
4
2
3
0
10
4 6 Iteration number
(a) η = 1 ×
8
10
2
4
×10
10−2
4 6 Iteration number
(b) η = 5 ×
8
10
2
4
×10
10−2
4 6 Iteration number
(c) η = 1 ×
8
10 4
×10
10−1
of
2
p ro
Fig. 4.3. Error histories of the ART, RK, PART and PRK methods for fancurvedtomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c). 9 8 7 6
ART RK PART PRK
ART RK PART PRK
8 7 6
4
Error
Error
5
Error
9
ART RK PART PRK
3
5 4
100 2
2
4 6 Iteration number
8
10
Pr e-
3
2
4
×10
(a) η = 1 × 10−2
4 6 Iteration number
8
10
2
4
×10
(b) η = 5 × 10−2
4 6 Iteration number
8
10 4
×10
(c) η = 1 × 10−1
Fig. 4.4. Error histories of the ART, RK, PART and PRK methods for sphericaltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
6
Error
Error
5 4
al
3
100
(a) η = 1 ×
8
10
urn
4 6 Iteration number
4
×10
10−2
2
7 6 5
3
4 6 Iteration number
(b) η = 5 ×
ART RK PART PRK
8
4
2
2
9
ART RK PART PRK
Error
9 8 7
ART RK PART PRK
10−2
8
10 4
×10
2
4 6 Iteration number
(c) η = 1 ×
8
10 4
×10
10−1
Fig. 4.5. Error histories of the ART, RK, PART and PRK methods for paralleltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
Jo
have an apparent gap in size. From these figures, we clearly see the semiconvergence of these algorithms and that the PRK method has a smaller minimum error than the ART, PART and RK methods. Then we track the errors kxk − x∗ k of the REK-ZF, REK-Du, REK-BW and PREK methods. In these methods, we need two initial guesses x0 and z0 , which are given to be zero vector and the right-hand side, respectively. The results are shown in Figure 4.6 – Figure 4.8 after 1.4 × 105 iterations. We can see that the behaviors of the PREK, REK-ZF, REK-Du and REK-BW methods break the semiconvergence horizon of the PRK method and the proposed algorithm outperforms the REK-ZF, REK-Du and REK-BW methods in terms of the iteration numbers. 5. Conclusions. For a constrained linear system, the PSIRT and PART methods are available, and the latter is observed to be faster than the former [17]. In this paper, we present the PRK method to further accelerate the convergence of PART, and develop its variant for the inconsistent case. We first give the PRK method to solve the constrained linear system with the noise-free right-hand side and analyze its convergence behavior. We show that the PRK method can converge exponentially in expectation to the solution for the consistent system, and find that the PRK method has a faster rate of convergence than the ART, PART
Journal Pre-proof 10
Projected randomized Kaczmarz method
1.5 REK-ZF REK-Du REK-BW PREK Error
Error
REK-ZF REK-Du REK-BW PREK
1
Error
0
10
REK-ZF REK-Du REK-BW PREK
1
0.5012
0.5 0.2512 10-1 0.1259 4
6 8 10 Iteration number
(a) η = 1 ×
12
2
14
4
4
×10
10−2
6 8 10 Iteration number
(b) η = 5 ×
12
14
2
4
4
×10
10−2
6 8 10 Iteration number
(c) η = 1 ×
of
2
12
14 4
×10
10−1
REK-ZF REK-Du REK-BW PREK
0
10
p ro
Fig. 4.6. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for fancurvedtomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c). 1.6 1.4 1.2
REK-ZF REK-Du REK-BW PREK
0
10
REK-ZF REK-Du REK-BW PREK
1
10-1
Error
Error
Error
0.8 0.6
0.4
-1
10 2
4
6 8 10 Iteration number
(a) η = 1 ×
12
14
Pr e-
0.2
-2
10
2
4
×10
10−2
4
6 8 10 Iteration number
(b) η = 5 ×
12
2
4
4
×10
10−2
6 8 10 Iteration number
(c) η = 1 ×
12
14 4
×10
10−1
Fig. 4.7. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for sphericaltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
REK-ZF REK-Du REK-BW PREK
100
1.6 1.4 1.2
REK-ZF REK-Du REK-BW PREK
1.2589
1 0.8 Error
al
Error
0.6310 Error
REK-ZF REK-Du REK-BW PREK
0.3162
0.6
0.4
0.1585
-1
1
2
3
urn
10
4 5 6 7 Iteration number
(a) η = 1 ×
8
10−2
9
10 4
×10
2
4
0.2 6 8 10 Iteration number
(b) η = 5 ×
10−2
12
14 4
×10
2
4
6 8 10 Iteration number
(c) η = 1 ×
12
14 4
×10
10−1
Fig. 4.8. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for paralleltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
Jo
and RK methods as exhibited in Figure 4.2. Then, we employ the PRK method to the constrained linear system with the noisy right-hand side (i.e., bn 6= 0). We observe that the PRK method will not converge to the LS solution for this case. As shown in Figure 4.3 - Figure 4.5, the iteration vector approaches a regularized solution in the early stage, while it is deteriorated by noise in the later iterations. Therefore, we present the semiconvergence analysis of the PRK method and give an upper estimate of the noise error. Finally, we obtain the PREK method by using the combination of the randomized orthogonal projection and PRK methods, which is a new LS solver for the inconsistent constrained linear system. We prove that the PREK method can converge exponentially in expectation to the LS solution for the inconsistent system. We can see that the PREK method outperforms the REK-ZF, REK-Du and REK-BW methods, which are revealed in Figure 4.6 – Figure 4.8. Moreover, like the successive over relaxation (SOR) method for the Gauss-Seidel (GS) method [20, 35], if we introduce an relaxation parameter λk to the PRK and PREK methods, we can have the following relaxed iterative schemes.
Journal Pre-proof 11
NIANCI WU AND HUA XIANG
³ ´ b −aT x (a) Relaxed PRK: xk = PC xk−1 + λk i kaii k2k−1 ai .
p ro
of
³ ´ Aj AT [b−zk ]i −aT i xk−1 (b) Relaxed PREK: zk = zk−1 − kAj kj2 zk−1 , xk = PC xk−1 + λk a . 2 i kai k Choosing appropriate relaxation parameters can improve the convergence of the relaxed PRK and PREK methods. We remark that the choice of probability distribution can significantly affect the performance of the PRK and PREK methods. In [18] a modified version of the RK method is presented to improve its convergence rate by utilizing the Johnson-Lindenstrauss dimension reduction technique. The other accelerated forms of the RK method are given in Bai and Wu’s series of work [2, 4]. They redefine an effective probability criterion for selecting the working rows from A and construct a greedy randomized Kaczmarz (GRK) method, which can strikingly improve the convergence rate of the original form, and its relaxed version is given in [4]. After that, another new greedy Kaczmarz method is also proposed in [44]. Based on these facts, we can employ the above strategies to the PRK and PREK methods and further accelerate their convergence. Acknowledgment. We thank the referees of an earlier version of this paper for several comments that lead to improvements of the theory as well as the presentation. REFERENCES
Jo
urn
al
Pr e-
[1] M.S. Andersen and P.C. Hansen. Generalized row-action methods for tomographic imaging. Numer. Alg., (2014)67: 121144. [2] Z. Z. Bai and W. T. Wu. On greedy randomized Kaczmarz method for solving large sparse linear systems. SIAM J. Sci. Comput., (2018)40: A592-A606. [3] Z. Z. Bai and W. T. Wu. On convergence rate of the randomized Kaczmarz method. Linear Algebra Appl., (2018)553: 252-269. [4] Z. Z. Bai and W. T. Wu. On relaxed greedy randomized Kaczmarz methods for solving large sparse linear systems. Appl. Math. Lett., (2018)83: 21-26. [5] Z. Z. Bai and W. T. Wu. On Partially randomized extended Kaczmarz method for solving large sparse overdetermined inconsistent linear systems. Linear Algebra Appl., (2019), https://doi.org/10.1016/j.laa.2019.05.005 [6] A. Baghban, et al., Developing an ANFIS based swarm concept model for estimating relative viscosity of nanofluids. Engineering Applications of Computational Fluid Mechanics, (2019)13: 26-39. [7] R. Bramley and A. Sameh. Row projection methods for large nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., (1992)13: 168-193. [8] Y. Censor, D. Gordon, and R. Gordon. Component averaging: An efficient iterative parallel algorithm for large and sparse unstructured problems. Parallel Comput., (2001)27: 777-808. [9] C.T. Cheng, et al.,Three-person multi-objective conflict decision in reservoir flood control. European Journal of Operational Research, (2002)142: 625-631. [10] G. Cimmino, Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari, La Ricerca Scientifica, XVI, Series II, Anno IX, (1938)1: 326-333. [11] J. W. Demmel. The probability that a numerical analysis problem is difficult. Math. Comp., (1988)50: 449-480. [12] K. Du. Tight upper bounds for the convergence of the randomized extended Kaczmarz and Gauss-Seidel algorithms. Numer. Linear Algebra Appl., (2019)26: e2233. [13] B. Eicke. Iteration methods for convexly constrained ill-posed problems in Hilbert space. Numer. Funct. Anal. Optim., (1992)13: 413-429. [14] T. Elfving, T. Nikazad, and P. C. Hansen. Semi-convergence and relaxation parameters for a class of SIRT algorithms. Electron. Trans. Numer. Anal., (2010)37: 321-336. [15] T. Elfving and P. C. Hansen. Unmatched projector back projector pairs perturbation and convergence analysis. SIAM. Sci. Comput., (2018)40: A573-A591. [16] T. Elfving, P. C. Hansen, and T. Nikazad. Semiconvergence and relaxation parameters for projected SIRT algorithms. SIAM J. Sci. Comput., (2012)34: A2000-A2017. [17] T. Elfving, T. Nikazad, and P. C. Hansen. Semi-convergence properties of Kaczmarz’s method. Inverse Problems, (2014)30: 055007. [18] Y. C. Eldar and D. Needell. Acceleration of randomized Kaczmarz method via the Johnson-Lindenstrauss Lemma. Numer. Alg., (2011)58: 163-177. [19] R. Escalante and M. Raydan. Alternating Projection Methods. SIAM Philadelphia, 2011. [20] G. H. Golub and C. V. Loan. Matrix Computations, 4th ed., Johns Hopkins University Press, Baltimore, MD, 2013. arik. Randomized iterative methods for linear systems. SIAM J. Matrix Anal. Appl., (2015)36: [21] R. M. Gower and P. Richt´ 1660-1690. [22] P. C. Hansen and M. S. Hansen. AIR Tools–A MATLAB package of algebraic iterative reconstruction methods. J. Comput. Appl. Math., (2012)236: 2167-2178. [23] P. C. Hansen and J. S. Jørgensen. AIR Tools II: algebraic iterative reconstruction methods, improved implementation. Numer. Alg. , (2018)79:107-137. [24] G. T. Herman. Fundamentals of Computerized Tomography: Image Reconstruction from Projections. 2nd edn (New York: Springer), 2009. [25] Y. Jiao, B. Jin, and X. Lu. Preasymptotic convergence of randomized Kaczmarz method. Inverse Problems, (2017)33: 125012. aherte Aufl¨ osung von Systemen linearer Gleichungen, in Bulletin International de l’Acad´ emie Polon[26] S. Kaczmarz. Angen¨ aise des Sciences et des Lettres. Classes des Sciences Mathmatiques et Naturelles. Series A: Sciences Mathmatiques, 1937: 355-357. [27] L. Landweber. An iterative formula for Fredholm integral equations of the first kind. Amer. J. Math., (1951)73: 615-624. [28] D. Leventhal and A. S. Lewis. Randomized methods for linear constraints: convergence rates and conditioning. Math.
Journal Pre-proof 12
Projected randomized Kaczmarz method
Jo
urn
al
Pr e-
p ro
of
Oper. Res., (2010)35: 641-654. [29] J. Lin and S. Li. Convergence of projected Landweber iteration for matrix rank minimization. Appl. Comput. Harm. Anal., (2014)36: 316-325. [30] R. Moazenzadeh, et al., Coupling a firefly algorithm with support vector regression to predict evaporation in northern Iran. Engineering Applications of Computational Fluid Mechanics, (2018)12: 584-597. [31] D. Needell. Randomized Kaczmarz solver for noisy linear systems. BIT Numer. Math., (2010)50: 395-403. [32] T. Nikazad, M. Abbasi, and T. Elfving. Error minimizing relaxation strategies in Landweber and Kaczmarz type iterations. J. Inverse Ill-Posed Probl., (2017)25: 35-56. [33] C. Popa. Least-squares solution of overdetermined inconsistent linear systems using Kaczmarz’s relaxation. Int. J. Comput. Math., (1995)55: 79-89. [34] C. Popa. Characterization of the solutions set of inconsistent least-squares problems by an extended Kaczmarz algorithm. Korean J. Comput. Appl. Math.,(1999)6: 51-64. [35] Y. Saad. Iterative Methods for Sparse Linear Systems. 2nd edition, SIAM, Philadelphia, PA, USA, 2000. [36] S. Samadianfard, et al., Daily global solar radiation modeling using data driven techniques and empirical equations in a semi-arid climate. Engineering Applications of Computational Fluid Mechanics, 2019 (13): 142-157. [37] R. W. Schafer, R. M. Mersereau, and M. A. Richards. Constrained iterative restoration algorithms. Proc. IEEE, (1981)69: 432-450. [38] T. Strohmer and R. Vershynin. A randomized Kaczmarz algorithm for linear systems with exponential convergence. J. Fourier Anal. Appl., (2009)15: 262-278. [39] C. Vogel. Computational Methods for Inverse Problems. SIAM, Philadelphia, 2002. [40] C.L. Wu, et al.,Rainfall-runoff modeling using artificial neural network coupled with singular spectrum analysis. J. Hydrol., ( 2011)399: 394-409. [41] H. Xiang and L. Zhang. Randomized iterative methods with alternating projections. arXiv preprint arXiv: 1708.09845, 2017. [42] Z.M. Yaseen, et al., An enhanced extreme learning machine model for river flow forecasting: state-of-the-art, practical applications in water resource engineering area and future research direction. J. Hydrol., (2019)569: 387-408. [43] H. Zhang and L. Z. Cheng. Projected Landweber iteration for matrix completion. J. Comput. Appl. Math., (2010)235: 593-601. [44] J. J. Zhang. A new greedy Kaczmarz algorithm for the solution of very large linear systems. Appl. Math. Letters, (2019)91: 207-212. [45] A. Zouzias and N. M. Freris. Randomized extended Kaczmarz for solving least-squares. SIAM J. Matrix Anal. Appl., (2012)34: 773-793. [46] X. Zheng, J. Yang, L. Li, et al. Modified projected landweber super-resolution algorithms for passive millimeter wave imaging. International Conference on Communications, 2008.
PII: DOI: Reference:
S0377-0427(19)30677-6 https://doi.org/10.1016/j.cam.2019.112672 CAM 112672
To appear in:
Journal of Computational and Applied Mathematics
Received date : 12 March 2019 Revised date : 10 August 2019 Please cite this article as: N. Wu and H. Xiang, Projected randomized Kaczmarz methods, Journal of Computational and Applied Mathematics (2020), doi: https://doi.org/10.1016/j.cam.2019.112672. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Elsevier B.V. All rights reserved.
Journal Pre-proof
PROJECTED RANDOMIZED KACZMARZ METHODS NIANCI WU
†
∗
AND HUA XIANG †
of
Abstract. For a constrained linear ill-posed problems, the projected algebra reconstruction technique (PART) usually has a faster convergence rate than the projected simultaneous iterative reconstruction technique (PSIRT), but can not be applied to the inconsistent case. In this work we propose the projected randomized Kaczmarz method (PRK) to further accelerate the convergence of PART and analyze its convergence. With the noise-free right-hand side, we show that the PRK method can converge exponentially in expectation to the solution for the consistent system. With the noisy right-hand side, we present some insights into the semiconvergence property of the PRK method and propose an extended version of PRK (PREK) method to seek the least squares solution. We prove that the PREK method can converge exponentially in expectation to the least squares solution for the inconsistent systems. Key words. projected algebra reconstruction technique, projected randomized iterative method, exponential convergence, semiconvergence
p ro
AMS subject classifications. 15A23, 65F30, 65F22, 65F10, 65R32
1. Introduction. Consider the constrained discrete ill-posed problems Ax = b subject to x ∈ C,
(1.1)
(1.2)
urn
al
Pr e-
where A ∈ Rm×n (m ≥ n), b ∈ Rm , and C is the nonempty closed convex set. The discrete ill-posed problems are characterized by the coefficient matrices with a very large condition number. This implies that the naive solution is very sensitive to any perturbation of the right-hand side, representing the errors in the data. The numerical solution of these systems requires the use of iterative regularization methods. Furthermore, iterative methods are usually more suitable than direct methods for a large scale problem. Some numerical techniques in the real-life applications can be found in different fields (see [6, 9, 30, 36, 40, 42]). It is often necessary to incorporate additional constrains to the iterative regularization methods, which can typically lead to smaller iteration errors in computational tomography as well as many other image problems (see [1, 14, 16, 17]). The common closed convex set includes nonnegativity constrain, box constrain, finite support constrain and bandlimiting constrain (see [37]). Assuming that S and T are two closed subspaces of a Hilbert space H, to find the best approximation to x from S ∩ T , the method of alternating projection (MAP) [19] suggests that we can first project x onto S, the obtained element is then projected onto T , and the algorithm projects alternatively onto S and T . In this way it generates a sequence of elements that converges to PS∩T (x). Combining the projection methods with the reconstruction techniques, e.g., simultaneous iterative reconstruction technique (SIRT) [14], in each step of the MAP method, Elfving, Hansen and Nikazad [16] give the following projected SIRT (PSIRT) method, which has extensive applications in image restoration, matrix rank minimization and matrix completion (see [16, 29, 39, 43, 46]). Algorithm PSIRT. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 21 = xk−1 + λk−1 AT M (b − Axk−1 ),
xk = PC (xk− 21 ),
Jo
where {λk }∞ k=0 is a sequence of positive relaxation parameters, M is a symmetric positive definite matrix and PC is the metric projection onto the nonempty closed convex set C. Some well-known methods can be written in the form (1.2) for appropriate choices of the matrix M . The following three choices of M are quite popular: (1) M = Im where Im is the identity Pn matrix of order m, 1 (2) M = m diag(1/kai k2 ) where aTi denotes the i-th row of A, and (3) M = diag(1/ j=1 Nj a2ij ) where Nj is the number of nonzeros in the j-th column of A and aij refers to the (i, j) entry of A. These three cases correspond to the classical projected Landweber method [27, 29, 43, 46], the projected Cimmino method [10, 16], and the projected component averaging (CAV) method [1, 8, 16], respectively. Let h be the spectral radius of the matrix AT M A and assume that 0 ≤ ² ≤ λk ≤ (2 − ²)/h. If ² > 0, then the iterates of SIRT with and without projections converge to a solution of min kAx − bkM [14], where kuk2M = uT M u is the norm induced by the symmetric positive definite matrix M for any vector u. Moreover, developing appropriate step size strategies can significantly improve the convergence rate of SIRT, for more details, refer to [14, 32]. The Kaczmarz method [26] proposed in 1937 is a simple yet powerful technique for solving the linear system. In tomography research it is also known as the algebra reconstruction technique (ART) [24]. ART ∗ This work was supported by the National Natural Science Foundation of China under Grant No. 11571265 and NSFC-RGC No. 11661161017. † School of Mathematics and Statistics, Wuhan University, Wuhan 430072, PR China. ([email protected]).
1
Journal Pre-proof 2
Projected randomized Kaczmarz method
is a particular case of row projection methods that suits for parallel computations and large scale problems since each step only requires one row of matrix A and no matrix-vector products [7], while SIRT makes 1 use of all the equations simultaneously. For example, by taking M = m diag(1/kai k2 ) and the relaxation parameter as constant 1, SIRT automatically reduces to the Cimmino method, which computes the next iterative by averaging all the projections of the previous iterates. The iterate schemes of the Cimmino and ART methods are given by
of
xk− 21
¶ m m µ 1 X 1 X bi − aTi xk−1 = PHi (xk−1 ) ≡ xk−1 + ai m i=1 m i=1 kai k2 ¶ m µ 1 X bi − aTi xk−1 a = xk−1 + AT M (b − Axk−1 ) = xk−1 + i m i=1 kai k2
and bi − aTi xk−1 ai kai k2
(i = 1, 2, · · · , m, 1, 2, · · · , m, 1, 2, · · · ),
p ro
xk− 12 = PHi (xk−1 ) ≡ xk−1 +
respectively, where PHi is the metric projection onto the hyperplane Hi = {x ∈ Rn : hai , xi = bi },
Pr e-
and bi is the i-th component of the right-hand side vector b ∈ Rm . In addition, we find that if we have B = AAT at the beginning, then ART only costs 4n + 1 flopping operations (flops), while the Cimmino method requires for 4mn + m flops, at each iteration step. In the recent work of Elfving, Hansen and Nikazad [17], they have found empirically that ART usually has a faster convergence rate than some SIRT such as the Landweber and Cimmino methods. The geometry of the Cimmino and ART methods are also shown in Figure 1.1 when m = 2, which can intuitively explain the fast convergence of the ART method compared to the Cimmino method. H1
xk-1
xk-1
al
H1
urn
H2
(a) The Cimmino method
H2
(b) The ART method
Fig. 1.1. The geometry of the Cimmino and ART methods. The hyperplane Hi = {x ∈ Rn : hai , xi = bi } for i = 1, 2.
(1.3)
Jo
By substituting SIRT with ART in the PSIRT method, Andersen and Hansen present the following projected ART (PART) method in [1]. Algorithm PART. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 12 = PHi (xk−1 ),
xk = PC (xk− 12 ),
i = 1, 2, · · · , m, 1, 2, · · · , m, 1, 2, · · · ,
where PHi and PC are the metric projections onto the hyperplane Hi and the nonempty closed convex set C, respectively. The theoretically convergence rate of ART is very challenging, which is reflected by the fact that the convergence rate of the method depends strongly on the ordering of the equations. Strohmer and Vershynin in an influential paper [38] use the rows of A in ART with a random order, called the randomized Kaczmarz (RK) method below, rather than the given order. This strategy can guarantee an exponential convergence rate1 of the RK method for consistent linear systems, where the convergence rate depends on the condition number. This result is then extended and refined in various directions [2, 3, 4, 5, 21, 28, 31, 41, 45]. 1 A randomized iterative method that computes a sequence of random vectors {x }∞ k k=0 is said to converge in expectation to a vector x∗ if and only if E[kxk − x∗ k2 ] → 0 as k → ∞, where the expectation is taken over the random choice of the methods.
Journal Pre-proof 3
NIANCI WU AND HUA XIANG
of
Inspired by the PART and RK methods, we propose the projected randomized Kaczmarz (PRK) method to solve the constrained linear ill-posed problems. First, the PRK method is given in section 2, including its convergence analyses for the noise-free and noisy right-hand sides. Second, we present an extended version of the PRK (PREK) method that can converge exponentially in expectation to the least square solution for the constrained inconsistent systems in section 3. Finally, we demonstrate the convergence behaviors of the proposed algorithms with the numerical examples from medical tomography. Throughout the paper k·k denotes the vector and matrix 2-norm, and k·kF denotes the matrix Frobenius norm. For a rectangular matrix A, we denote by A† the Moore-Penrose pseudoinverse of A [20], R(A) the range of A and N (AT ) the null space of AT . Given any vector b ∈ Rm , we can uniquely write it as b = br + bn , where br = AA† b and bn = (Im − AA† )b are the projections of b onto R(A) and N (AT ), respectively. We indicate by Ek−1 [·] the expected conditional on the k-th iteration, and from the law of the iterated expectation we have E[Ek−1 [·]] = E[·]. Contributions. We comment on the main contributions of this work as follows.
Pr e-
p ro
(i) A new scheme. The PART is a popular way to solve the constrained linear system, which is under the framework of the MAP method. This method suggests that at the first-half step, we use ART to the linear system Ax = b, and at the second-half step, the obtained element is projected onto the set C. Then this process continues alternatively. Based on the fact that the RK method has a faster convergence rate than ART, we replace ART with RK in the PART method and get the PRK method to solve the constrained linear system. We show that the PRK method possesses the expected exponential convergence rate. (ii) Semiconvergence analysis. We apply the PRK method to the case where the constrained linear system is inconsistent, which has the noisy right-hand side (i.e., bn 6= 0). We observe that the PRK method will not converge to the least squares (LS) solution of the inconsistent system. Some insights into its semiconvergence property are presented also. To investigate how the noisy part bn is propagated during the PRK iterations, we split the total error into iteration error and noise error, and give an upper estimate of the noise error. (iii) A extended variant for the LS solution. We utilize the randomized orthogonal projection method to find the approximation of bn and employ the PRK method to the approximate consistent system Ax ≈ b−bn . Then the PREK method is obtained, which is a new LS solver of the inconsistent constrained linear system. We prove that the PREK method can converge exponentially in expectation to the LS solution for the inconsistent system. 2. The projected randomized Kaczmarz method. In this section, we first present the PRK method to solve the constrained consistent linear ill-posed problems. Then, its convergence analyses are followed for the noise-free and noisy right-hand sides separately.
(2.1)
urn
al
2.1. The PRK method. Motivated by the PART process, the PRK method is mainly based on the next two generic step procedures. We first project the iterate vector xk−1 onto the i-th hyperplane Hi in a random order, and then project the current vector onto the closed convex set C. Continuing with this process alternatively, the PRK method can be formulated as follows. Algorithm PRK. Given an initial guess x0 ∈ Rn , for k = 1, 2, · · · until {xk }∞ k=0 converges, compute xk− 21 = PHi (xk−1 ),
xk = PC (xk− 12 ),
Jo
where PHi and PC are the metric projections onto the hyperplane Hi and the nonempty closed convex set C, respectively, the row index i is chosen from the set [m] := {1, 2, · · · , m} at random with probability pi = kai k2 /kAk2F . Remark 1. Let Sk = {x ∈ Rn : Aτk x = bτk , τk ⊂ [m]}, where bτk ∈ Rτk , Aτk ∈ Rτk ×n is formed by the rows of A, then for any y ∈ / Sk the projection onto Sk is given by PSk (y) = y + ATτk (Aτk ATτk )† (bτk − Aτk y).
Instead of just using one row aTi of A in (2.1), we can work on several rows simultaneously, it is natural to extend the PRK method to block variants. Different from ART, which selects the row index cyclically in a deterministic fashion, the RK method selects the row index at random with a probability criterion. Thus, the PRK method can be viewed as the random version of the PART method. Additionally, at the first-half iteration step, the PART and PRK methods need the same flops. From the geometry viewpoint, we know that if we project the current vector xk− 12 onto the hyperplane Hi (i = 1, 2, · · · , m) as far as possible at the first-half step of the PART method, then it could lead to a faster convergence rate. This also can be seen from the Figure 1.1(b). In designing a random version of the PART method, i.e., the PRK method, we aim at grasping larger entries of the residual vector and selecting the working row according to a probability criterion, rather than in their given order. It is possible to avoid a hyperplane, which has a small distance from xk−1 , to be selected.
Journal Pre-proof 4
Projected randomized Kaczmarz method
k I 2 E[keIk k2 ] ≤ (1 − κ−2 A ) · ke0 k ,
p ro
(2.2)
of
2.2. Convergence analysis of PRK with noise-free right-hand side. Let x∗ be the solution of (1.1) and xk denote the k-th iterate of Algorithm PRK using the noise-free right-hand side, i.e., b = Ax∗ , then bi = aTi x∗ for i = 1, 2, · · · , m. Define the iteration error eIk := xk − x∗ , which is also known as data error, see [14, 16, 17]. Next, we will give an upper estimate on E[keIk k2 ] for the PRK method. This represents our first main theoretical result. Lemma 2.1. [19, Proposition 4.1] If C is a closed and convex set in any Hilbert space H, then the projection operator PC : H → C is non-expansive, i.e., for any elements α, β ∈ H, kPC (α) − PC (β)k ≤ kα − βk. By utilizing the above lemma, we can deduce the convergence property of the PRK method with noisefree right-hand side and establish the following theorem. Theorem 2.2. Let the constrained linear system (1.1), with the coefficient matrix A ∈ Rm×n and the right-hand side b ∈ Rm , be consistent. Then from any initial guess x0 ∈ R(AT ), the iteration sequence ∗ {xk }∞ k=0 generated by the PRK method converges to the solution x in expectation. Moreover, the mean squared iteration error for k = 0, 1, 2, · · · satisfies
where the generalized condition number κA = kAkF · kA† k. Proof. Based on the non-expansivity of the projection operator in Lemma 2.1, it holds that m m X X kai k2 kai k2 ∗ 2 ∗ 2 = · kPC [PHi (xk−1 )] − PC (x )k ≤ 2 2 · kPHi (xk−1 ) − x k kAk kAk F F i=1 i=1
= (2.3) =
m X kai k2 bi − aTi xk−1 ∗ · k(x − x ) + ai k2 k−1 2 2 kAk ka k i F i=1
Pr e-
Ek−1 [keIk k2 ]
m X kai k2 ai aTi · k(In − )eIk−1 k2 2 2 kAk ka k i F i=1
m X kai k2 ai aTi AAT I I I T · (eIk−1 )T (In − )e = )e . = (e ) (I − n k−1 k−1 2 kAkF kai k2 kAk2F k−1 i=1
al
By induction, xk−1 is in the row space of A for any k when x0 ∈ R(AT ). In addition, the same is true for x∗ by the definition of psedoinverse. This means that eIk−1 = xk−1 − x∗ is also in the row space of A. Therefore, we can have µ ¶ σr2 (A) AAT I I 2 I T )e ≤ 1− · keIk−1 k2 = (1 − κ−2 (ek−1 ) (In − A ) · kek−1 k . kAk2F k−1 kAk2F
urn
By taking the full expectation for both sides of the inequality (2.3), we obtain I 2 E[keIk k2 ] ≤ (1 − κ−2 A ) · E[kek−1 k ].
Jo
Repeating the same argument k − 1 times we get the estimate (2.2). Strohmer and Vershynin in [38] prove a similar exponential upper bound on the expected convergence rate for the noise-free right-hand side about the RK method. Theorem 2.2 indicates that the convergence rate of the PRK method is not less than that of the RK method due to the non-expansivity of the projection operator. Besides, the expected convergence rate of the RK method in [38] is 1−b κ−2 , where κ b = kAkF ·kA−1 k is the scaled condition number introduced by Demmel [11]. It is based on the restrictive assumption that A has full column rank, while the convergence rate of the PRK method only depends on the generalized condition number κA without imposing this assumption. In practice, one usually observes the typical semiconvergence phenomenon for iterative regularization methods, e.g., unprojected/projected SIRT in [16, 14], ART in [1, 17] and the general iteration scheme in [15]. That is, the iterative procedure of these methods are damaged by noise. The iteration vector approaches a regularized solution in the early stage, while it is deteriorated by noise in the later iterations. Next we will study the semiconvergence behavior of the PRK method. 2.3. Semiconvergence analysis of PRK with noisy right-hand side. If the constrained linear system (1.1) is inconsistent, we use x∗ = A† b ∈ C to represent its LS solution. Then it holds that the ek denote the k-th iterates of the PRK method using the noisy right-hand side b = Ax∗ + bn . Let xk and x ek − x∗ and the noise error noise-free and noisy right-hand sides, respectively. Define the total error e ek := x N ek − xk , where the noise error is also known as approximation error, see [14, 16, 17]. ek := x
Journal Pre-proof 5
NIANCI WU AND HUA XIANG
I Clearly, the total error e ek satisfies that e ek = eN k + ek . Classical convergence theorem usually focuses on eIk (e.g., Theorem 2.2), while assuming that eN is negligible. However, this is not the case for ill-posed and k noisy problems. The interplay of these two error terms determines the final total error. We next investigate how the noisy part bn in the right-hand side is propagated during the PRK iterations. Theorem 2.3. Let the linear system (1.1) be with the coefficient matrix A ∈ Rm×n . The iteration sequences {xk }∞ xk }∞ k=0 and {e k=0 are generated by the PRK method with noise-free right-hand side b = m br ∈ R and noisy right-hand side b = br + bn ∈ Rm , respectively. Assuming that the initial guesses e0 = x0 ∈ R(AT ), then the mean squared noise error eN ek − xk is bounded by x k =x
√
k kbn k2 · , σr (A) kAkF
where σr (A) is the smallest singular value of A Proof. Direct calculation yields that µ
ek−1 + x
ek−1 eTi (br + bn ) − aTi x ai 2 kai k
¶
µ
¶ eTi bn − aTi xk−1 a . i kai k2
p ro
eN k = PC
of
2 E[keN k k ]≤
(2.4)
− PC
xk−1 +
Using again the non-expansivity of the projection operator, we can have
m m X X kai k2 kai k2 ai aTi eTi bn N 2 N · ke · k(I − k ≤ )e + a k2 n k k−1 2 2 2 2 i kAk kAk ka k ka k i i F F i=1 i=1 · ¸ m 2 X ai aTi (eTi bn )2 kai k N T N · (e ) (I − )e + = n k−1 kAk2F kai k2 k−1 kai k2 i=1 "m # m X kai k2 X ai aTi (eTi bn )2 T N = (eN · (I − ) e + n k−1 ) k−1 kAk2F kai k2 kAk2F i=1 i=1
Pr e-
2 Ek−1 [keN k k ]=
N 2 ≤ (1 − κ−2 A ) · kek−1 k +
kbn k2 , kAk2F
where ei is the i-th column of identity matrix Im with order m. By taking the full expectation for both sides of the above inequality, it indicates that k−1 2 X kbn k2 −2 k −2 j kbn k N 2 ≤ (1 − κ ) · ke k + . (1 − κ ) · 0 A A kAk2F kAk2F j=0
al
−2 2 N 2 E[keN k k ] ≤ (1 − κA ) · E[kek−1 k ] +
urn
Assuming x = 1/κA ∈ (0, 1), the above inequality reads 2 E[keN k k ]≤
k−1 X j=0
j (1 − κ−2 A ) ·
e0 = x0 . Let y = 1 − x2 ∈ (0, 1), then since x
Jo
(2.5)
1 − (1 − x2 )k = x
s
(1 − y k )2 ≤ 1−y
kbn k2 1 − (1 − x2 )k kbn k2 = · κ · , A kAk2F x kAk2F
s
p √ 1 − yk = 1 + y + · · · + y k−1 ≤ k, 1−y
and the bound (2.4) comes out. This completes the proof. √ Note that the bound in Theorem 2.3 includes a factor k, which is similar to the formula for the propagation of the noise error in ART and SIRT, see [14, 17], and is valid also for their projected versions, see [16, 17]. The iteration number k and the smallest singular value σr (A) of A play the role of the regularized parameters. The bound (2.4) might be quite pessimistic since it holds for any closed convex set C and the estimates used in (2.5) are rather crude. As stated in Theorem 2.3, when the system (1.1) is perturbed by noise or no longer consistent, it is a pity that the PRK method will not converge to the LS solution of the inconsistent system. It is because that at the k-th iterate the PRK method projects the current vector xk−1 onto the hyperplane aTi x = eTi b = eTi (br +bn ) rather than aTi x = eTi br for i = 1, · · · , m. To resolve this problem, we provide the following variant of the PRK method.
Journal Pre-proof 6
Projected randomized Kaczmarz method
zk = zk−1 −
(3.1)
Aj ATj zk−1 , kAj k2
of
3. The projected randomized extended Kaczmarz method. In this section we will present a projected randomized iterative LS solver that converges in expectation to the LS solution of the inconsistent linear system with constraints. For the general inconsistent linear system, Zouzias and Freris in [45] modify the RK method and obtain a random version of extended Kaczmarz method motivated by the work of Popa in [33, 34]. We call it REKZF. One step of the REK-ZF method consists of two stages. The first stage uses one step RK update for the linear system Ax = b − zk−1 and generates xk . The second stage applies one step randomized orthogonal projection update to the homogeneous linear system AT z = 0 and outputs zk , which admits the following iterate scheme
Pr e-
p ro
where Aj is the j-th column of A and the index j is chosen at random. Actually, the iterate vector zk is the projection of zk−1 onto the hyperplane ATj z = 0, and it is equivalent to apply one RK update to AT z = 0. The REK-ZF method has an expected exponential convergence rate and immediately results in many works devoted to various aspects of the new algorithm, which includes distinct bounds on their performance and its extensions. Du presents a slightly different REK-ZF (denoted by REK-Du) method in [12], which generates xk by one-step RK update for the linear system Ax = b − zk (used in REK-Du) instead of Ax = b − zk−1 (used in REK-ZF) from xk−1 , because zk is a better approximation of bn than zk−1 . Recently, Bai and Wu [5] compute zk from (3.1) in a given cyclic order rather than in a dynamic random order, which can obtain the approximation of bn more accurately, and present the partial REK-ZF method. We call it REK-BW. It is proved that the REK-Du and REK-BW methods all have the smaller upper bounds for the expected convergence rate than that of the REK-ZF method (see [12, Thm.2] and [5, Thm.3.1]). 3.1. The PREK method. We know from Theorem 2.3 that the reason why the PRK method fails to converge to LS solution of the inconsistent problems is that bn 6= 0. Then, we split the PREK method into two computational stages. The first stage outputs zk by utilizing one step randomized orthogonal projection update to the linear system AT z = 0. The second stage uses one step PRK update to the approximate consistent linear system Ax = b − zk with x ∈ C. Note that if the nonempty closed convex set C = Rn , the PREK method could reduce to the REK-Du method. The pseudo-code of the PREK method is given in the following algorithm. Algorithm PREK. Given two initial guesses x0 ∈ Rm and z0 ∈ b + R(A), for k = 1, 2, · · · until ∞ {xk }k=0 converges, compute zk = PHe j (zk−1 ),
xk− 12 = PHb i (xk−1 ),
al
(3.2)
xk = PC (xk− 21 ),
where PHe j , PHb i and PC are the metric projections onto the hyperplanes
urn
e j = {z ∈ Rm : hAj , zi = 0}, H
b i = {x ∈ Rn : hai , xi = eT (b − zk )} H i
and the nonempty closed convex set C, respectively, the column index j is chosen from the set [n] at random with probability qj = kAj k2 /kAk2F and the row index i is chosen from the set [m] at random with probability pi = kai k2 /kAk2F .
(3.3)
Jo
3.2. Convergence analysis of PREK. For the convergence property of the PREK method, we can establish the following theorem. First, we review a known result needed in the remaining parts of this section. Lemma 3.1. [12, Thm.1] Let A ∈ Rm×n and b ∈ Rm . Let zk denote the k-th iterate of the randomized orthogonal projection method applied to AT z = 0 with initial guess z0 ∈ b + R(A). The iteration sequence {zk }∞ k=0 converges to the solution br in expectation. Moreover, the mean squared error for k = 0, 1, 2, · · · satisfies E[kzk − bn k2 ] ≤ ρk kz0 − bn k2 ,
where ρ = 1 − σr2 (A)/kAk2F = 1 − κ−2 A . Theorem 3.2. Let the constrained linear system (1.1), with the coefficient matrix A ∈ Rm×n and the right-hand side b = br + bn ∈ Rm , be inconsistent. Then from any initial guesses x0 ∈ R(AT ) and z0 ∈ b + R(A), the iteration sequence {xk }∞ k=0 generated by the PREK method converges to the solution x∗ = A† b in expectation. Moreover, the mean squared error for k = 0, 1, 2, · · · satisfies (3.4)
E[kxk − x∗ k2 ] ≤
kρk kz0 − bn k2 + ρk kx0 − x∗ k2 . kAk2F
Journal Pre-proof 7
NIANCI WU AND HUA XIANG
Proof. After the k-th iteration of the PREK method, it yields that kxk − x∗ k2 = kPC (xk−1/2 ) − PC (x∗ )k2
eTi (b − zk ) − aTi xk−1 ai − A† bk2 kai k2 eT (bn − zk ) eT br − aTi xk−1 = kxk−1 + i ai − A† b + i ai k2 2 kai k kai k2 eT AA† b − aTi xk−1 eT (bn − zk ) = kxk−1 + i ai − A† b + i ai k2 2 kai k kai k2 ¶ µ eT (bn − zk ) ai aTi (xk−1 − A† b) + i = k In − a i k2 2 kai k kai k2 ¶ µ eTi (bn − zk ) ai aTi † 2 = k In − (x − A b)k + k a i k2 , k−1 kai k2 kai k2 ≤ kxk−1/2 − x∗ k2 = kxk−1 +
p ro
of
(3.5)
(3.6)
Pr e-
where the second line comes from the the non-expansivity of the projection operator, the´third line follows ³ ai aT T from the fact that b = bn + br and the last line is from the orthogonality ai In − kai ki2 = 0. Like [12], we use Ek−1 [·] to denote the expected value conditional on the first k−1 iterations of the PREK method, i.e., Ek−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 ], where j` and i` (` = 1, · · · , k−1) mean that the `-th column and row are chosen as Ejk−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 , ik ] and Eik−1 [·] = E[·|j1 , i1 , · · · , jk−1 , ik−1 , jk ], respectively. Then, by the law of total expectation, i.e., Ek−1 [·] = Ejk−1 [Eik−1 [·]], it follows that ¸ ¸¸ · · ¸¸ · T · T · T e (bn − zk ) (ei (bn − zk ))2 ei (bn − zk ) j j 2 2 i i a k = E a k = E E E k Ek−1 k i i i k−1 k−1 k−1 k−1 kai k2 kai k2 kai k2 # "m · ¸ X kai k2 (eT (bn − zk ))2 kbn − zk k2 j j i · = Ek−1 = Ek−1 kAk2F kai k2 kAk2F i=1 =
£ ¤ 1 ρk 2 E kb − z k ≤ kbn − z0 k2 , k−1 n k kAk2F kAk2F
al
where the last line is from Lemma 3.1. Similar to the derivation of Theorem 2.2, we have that (3.7) · µ ¶ ¸ £ ¤ £ ¤ ai aTi † 2 † 2 (xk−1 − A b)k ≤ (1 − κ−2 = ρ Ek−1 kxk−1 − A† bk2 . Ek−1 k In − A ) · Ek−1 kxk−1 − A bk 2 kai k
urn
In accordance with formulas (3.6) and (3.7), the formula (3.5) then further results in the estimate · T ¸ · µ ¶ ¸ £ ¤ e (bn − zk ) ai aTi 2 † 2 Ek−1 kxk − x∗ k2 ≤ Ek−1 k i a k + E k I − (x − A b)k i k−1 n k−1 kai k2 kai k2 £ ¤ ρk ≤ kbn − z0 k2 + ρ Ek−1 kxk−1 − A† bk2 . 2 kAkF By taking the full expectation on both sides of the above inequality, it follows that
Jo
£ ¤ E kxk − x∗ k2 ≤
£ ¤ ρk kbn − z0 k2 + ρ E kxk−1 − A† bk2 . kAk2F
Then we immediately obtain the estimate in (3.4) by induction on the iteration index k. 4. Numerical experiments. In this section we investigate some numerical tests taken from the field of tomography image reconstruction. To discretize the tomography problems, we should divide the square domain into a grid of N × N pixels numbered from 1 to n (n = N 2 ). The j-th cell is assigned a constant value xj of the attenuation or travel time for j = 1, 2, · · · , n, such that the unknown vector x is a discrete version of the sought function, i.e., x = [x1 , · · · , xn ]T . Define the coefficient matrix A = [aij ] ∈ Rm×n , where aij is the length of the i-th ray through the j-th cell. Let the number of rays be m, then the measurement bi of the attenuation or travel time for the i-th ray can be written as bi =
n X j=1
aij xj = aTi x for i = 1, 2, · · · , m.
Journal Pre-proof 8
Projected randomized Kaczmarz method
For more details, refer to [22]. As the set C we use the positive orthant such that we enforce nonnegativity on the constructions. The orthogonal projector PC onto the positive orthant can be defined as follows. ½ yi , yi ≥ 0, [PC (y)]i = 0, yi < 0,
Pr e-
p ro
of
for any vector y = [y1 , · · · , yn ]T . Thus, a constrained discrete ill-posed problem is determined. Consider the Sheep-Logan phantom, its exact solution is shown in Figure 4.1(a). Based on various scanning modes, e.g., the fan-beam, spherical radon and parallel-beam method, where their geometries are shown with a few of rays in Figure 4.1((b)-(d)), the MATLAB package (AIR Tools, available from https://github.com/jakobsj/AIRToolsII) in [22, 23] presents three 2-dimension (2D) tomography test examples. They are fancurvedtomo, sphericaltomo and paralleltomo. These examples generate a sparse matrix A, an exact solution x∗ and an exact right-hand side b = Ax∗ . It gives us many flexibilities to adjust the input parameters, such as the size of the discrete domain (N ), the number of rays (p), the projection angles in degrees (θ) and the number of sources (s), where m = ps and n = N 2 . Therefore, we can assign values to these parameters and obtain various specific instantiations about the constrained linear systems. In this work, we typically set N = 40, p = 40, θ = 0 : 2 : 358 and s = 180. The resulting matrix is of size 7200 × 1600. All computations are started from the initial guess x0 = 0, and terminated once the relative solution error at the k-th iterate satisfies kxk − x∗ k ≤ kx∗ k · 10−6 or the number of iteration steps exceeds 105 . The numerical tests are performed on a Founder desktop PC with Intel(R) Core(TM) i5-7500 CPU 3.40 GHz by MATLAB R2016(a) with a machine precision of 10−16 .
(a) Sheep-Logan Phantom
(b) fan-beam
(c) spherical radon
(d) parallel-beam
al
Fig. 4.1. The exact solution of the Sheep-Logan phantom (a). The geometries of fan-beam (b), spherical radon (c) and parallel-beam (d), with a few of rays.
1
1
10
10
ART RK PART PRK Error
Error
100
ART RK PART PRK
100
urn
Error
ART RK PART PRK
-1
10
10-2
100
10-1
-3
10
2
4 6 Iteration number
Jo
(a) fancurvedtomo
8
10
4
×10
2
4 6 Iteration number
(b) sphericaltomo
8
10 4
×10
1
2
3
4 5 6 7 Iteration number
8
9
10 4
×10
(c) paralleltomo
Fig. 4.2. Error histories of the ART, RK, PART and PRK methods for fancurvedtomo(a), sphericaltomo (b) and paralleltomo (c) with noise-free right-hand side.
Case 1: noise-free right-hand side. In this case, we utilize the ART, RK, PART and PRK methods to solve the constrained consistent linear system. The errors kxk − x∗ k versus iteration number k for both methods with noise-free right-hand side are shown in Figure 4.2. In this figure, we can observe that the PRK method converges faster than the ART, PART and RK methods. Case 2: noisy right-hand side. To examine the semiconvergence phenomenon of the PRK method, like [17, 25], the noisy right-hand side b is defined as b = A · x∗ + η · kbk · ξ/kξk, where η is the relative noise level and ξ is a random vector generated by the MATLAB function randn. In this case the system (1.1) is no longer consistent, we first depict the convergence behavior of the ART, RK, PART and PRK methods in Figure 4.3 – Figure 4.5 with three different relative noise levels η = 1 × 10−2 , 5 × 10−2 and 1 × 10−1 , which
Journal Pre-proof 9
NIANCI WU AND HUA XIANG
9 8 7
ART RK PART PRK
9
ART RK PART PRK
6
7 6 Error
Error
5 Error
ART RK PART PRK
8
4
5
3
4
2
3
0
10
4 6 Iteration number
(a) η = 1 ×
8
10
2
4
×10
10−2
4 6 Iteration number
(b) η = 5 ×
8
10
2
4
×10
10−2
4 6 Iteration number
(c) η = 1 ×
8
10 4
×10
10−1
of
2
p ro
Fig. 4.3. Error histories of the ART, RK, PART and PRK methods for fancurvedtomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c). 9 8 7 6
ART RK PART PRK
ART RK PART PRK
8 7 6
4
Error
Error
5
Error
9
ART RK PART PRK
3
5 4
100 2
2
4 6 Iteration number
8
10
Pr e-
3
2
4
×10
(a) η = 1 × 10−2
4 6 Iteration number
8
10
2
4
×10
(b) η = 5 × 10−2
4 6 Iteration number
8
10 4
×10
(c) η = 1 × 10−1
Fig. 4.4. Error histories of the ART, RK, PART and PRK methods for sphericaltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
6
Error
Error
5 4
al
3
100
(a) η = 1 ×
8
10
urn
4 6 Iteration number
4
×10
10−2
2
7 6 5
3
4 6 Iteration number
(b) η = 5 ×
ART RK PART PRK
8
4
2
2
9
ART RK PART PRK
Error
9 8 7
ART RK PART PRK
10−2
8
10 4
×10
2
4 6 Iteration number
(c) η = 1 ×
8
10 4
×10
10−1
Fig. 4.5. Error histories of the ART, RK, PART and PRK methods for paralleltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
Jo
have an apparent gap in size. From these figures, we clearly see the semiconvergence of these algorithms and that the PRK method has a smaller minimum error than the ART, PART and RK methods. Then we track the errors kxk − x∗ k of the REK-ZF, REK-Du, REK-BW and PREK methods. In these methods, we need two initial guesses x0 and z0 , which are given to be zero vector and the right-hand side, respectively. The results are shown in Figure 4.6 – Figure 4.8 after 1.4 × 105 iterations. We can see that the behaviors of the PREK, REK-ZF, REK-Du and REK-BW methods break the semiconvergence horizon of the PRK method and the proposed algorithm outperforms the REK-ZF, REK-Du and REK-BW methods in terms of the iteration numbers. 5. Conclusions. For a constrained linear system, the PSIRT and PART methods are available, and the latter is observed to be faster than the former [17]. In this paper, we present the PRK method to further accelerate the convergence of PART, and develop its variant for the inconsistent case. We first give the PRK method to solve the constrained linear system with the noise-free right-hand side and analyze its convergence behavior. We show that the PRK method can converge exponentially in expectation to the solution for the consistent system, and find that the PRK method has a faster rate of convergence than the ART, PART
Journal Pre-proof 10
Projected randomized Kaczmarz method
1.5 REK-ZF REK-Du REK-BW PREK Error
Error
REK-ZF REK-Du REK-BW PREK
1
Error
0
10
REK-ZF REK-Du REK-BW PREK
1
0.5012
0.5 0.2512 10-1 0.1259 4
6 8 10 Iteration number
(a) η = 1 ×
12
2
14
4
4
×10
10−2
6 8 10 Iteration number
(b) η = 5 ×
12
14
2
4
4
×10
10−2
6 8 10 Iteration number
(c) η = 1 ×
of
2
12
14 4
×10
10−1
REK-ZF REK-Du REK-BW PREK
0
10
p ro
Fig. 4.6. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for fancurvedtomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c). 1.6 1.4 1.2
REK-ZF REK-Du REK-BW PREK
0
10
REK-ZF REK-Du REK-BW PREK
1
10-1
Error
Error
Error
0.8 0.6
0.4
-1
10 2
4
6 8 10 Iteration number
(a) η = 1 ×
12
14
Pr e-
0.2
-2
10
2
4
×10
10−2
4
6 8 10 Iteration number
(b) η = 5 ×
12
2
4
4
×10
10−2
6 8 10 Iteration number
(c) η = 1 ×
12
14 4
×10
10−1
Fig. 4.7. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for sphericaltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
REK-ZF REK-Du REK-BW PREK
100
1.6 1.4 1.2
REK-ZF REK-Du REK-BW PREK
1.2589
1 0.8 Error
al
Error
0.6310 Error
REK-ZF REK-Du REK-BW PREK
0.3162
0.6
0.4
0.1585
-1
1
2
3
urn
10
4 5 6 7 Iteration number
(a) η = 1 ×
8
10−2
9
10 4
×10
2
4
0.2 6 8 10 Iteration number
(b) η = 5 ×
10−2
12
14 4
×10
2
4
6 8 10 Iteration number
(c) η = 1 ×
12
14 4
×10
10−1
Fig. 4.8. Error histories of the REK-ZF, REK-Du, REK-BW and PREK methods for paralleltomo with three different relative noise levels η = 1 × 10−2 (a), 5 × 10−2 (b) and 1 × 10−1 (c).
Jo
and RK methods as exhibited in Figure 4.2. Then, we employ the PRK method to the constrained linear system with the noisy right-hand side (i.e., bn 6= 0). We observe that the PRK method will not converge to the LS solution for this case. As shown in Figure 4.3 - Figure 4.5, the iteration vector approaches a regularized solution in the early stage, while it is deteriorated by noise in the later iterations. Therefore, we present the semiconvergence analysis of the PRK method and give an upper estimate of the noise error. Finally, we obtain the PREK method by using the combination of the randomized orthogonal projection and PRK methods, which is a new LS solver for the inconsistent constrained linear system. We prove that the PREK method can converge exponentially in expectation to the LS solution for the inconsistent system. We can see that the PREK method outperforms the REK-ZF, REK-Du and REK-BW methods, which are revealed in Figure 4.6 – Figure 4.8. Moreover, like the successive over relaxation (SOR) method for the Gauss-Seidel (GS) method [20, 35], if we introduce an relaxation parameter λk to the PRK and PREK methods, we can have the following relaxed iterative schemes.
Journal Pre-proof 11
NIANCI WU AND HUA XIANG
³ ´ b −aT x (a) Relaxed PRK: xk = PC xk−1 + λk i kaii k2k−1 ai .
p ro
of
³ ´ Aj AT [b−zk ]i −aT i xk−1 (b) Relaxed PREK: zk = zk−1 − kAj kj2 zk−1 , xk = PC xk−1 + λk a . 2 i kai k Choosing appropriate relaxation parameters can improve the convergence of the relaxed PRK and PREK methods. We remark that the choice of probability distribution can significantly affect the performance of the PRK and PREK methods. In [18] a modified version of the RK method is presented to improve its convergence rate by utilizing the Johnson-Lindenstrauss dimension reduction technique. The other accelerated forms of the RK method are given in Bai and Wu’s series of work [2, 4]. They redefine an effective probability criterion for selecting the working rows from A and construct a greedy randomized Kaczmarz (GRK) method, which can strikingly improve the convergence rate of the original form, and its relaxed version is given in [4]. After that, another new greedy Kaczmarz method is also proposed in [44]. Based on these facts, we can employ the above strategies to the PRK and PREK methods and further accelerate their convergence. Acknowledgment. We thank the referees of an earlier version of this paper for several comments that lead to improvements of the theory as well as the presentation. REFERENCES
Jo
urn
al
Pr e-
[1] M.S. Andersen and P.C. Hansen. Generalized row-action methods for tomographic imaging. Numer. Alg., (2014)67: 121144. [2] Z. Z. Bai and W. T. Wu. On greedy randomized Kaczmarz method for solving large sparse linear systems. SIAM J. Sci. Comput., (2018)40: A592-A606. [3] Z. Z. Bai and W. T. Wu. On convergence rate of the randomized Kaczmarz method. Linear Algebra Appl., (2018)553: 252-269. [4] Z. Z. Bai and W. T. Wu. On relaxed greedy randomized Kaczmarz methods for solving large sparse linear systems. Appl. Math. Lett., (2018)83: 21-26. [5] Z. Z. Bai and W. T. Wu. On Partially randomized extended Kaczmarz method for solving large sparse overdetermined inconsistent linear systems. Linear Algebra Appl., (2019), https://doi.org/10.1016/j.laa.2019.05.005 [6] A. Baghban, et al., Developing an ANFIS based swarm concept model for estimating relative viscosity of nanofluids. Engineering Applications of Computational Fluid Mechanics, (2019)13: 26-39. [7] R. Bramley and A. Sameh. Row projection methods for large nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., (1992)13: 168-193. [8] Y. Censor, D. Gordon, and R. Gordon. Component averaging: An efficient iterative parallel algorithm for large and sparse unstructured problems. Parallel Comput., (2001)27: 777-808. [9] C.T. Cheng, et al.,Three-person multi-objective conflict decision in reservoir flood control. European Journal of Operational Research, (2002)142: 625-631. [10] G. Cimmino, Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari, La Ricerca Scientifica, XVI, Series II, Anno IX, (1938)1: 326-333. [11] J. W. Demmel. The probability that a numerical analysis problem is difficult. Math. Comp., (1988)50: 449-480. [12] K. Du. Tight upper bounds for the convergence of the randomized extended Kaczmarz and Gauss-Seidel algorithms. Numer. Linear Algebra Appl., (2019)26: e2233. [13] B. Eicke. Iteration methods for convexly constrained ill-posed problems in Hilbert space. Numer. Funct. Anal. Optim., (1992)13: 413-429. [14] T. Elfving, T. Nikazad, and P. C. Hansen. Semi-convergence and relaxation parameters for a class of SIRT algorithms. Electron. Trans. Numer. Anal., (2010)37: 321-336. [15] T. Elfving and P. C. Hansen. Unmatched projector back projector pairs perturbation and convergence analysis. SIAM. Sci. Comput., (2018)40: A573-A591. [16] T. Elfving, P. C. Hansen, and T. Nikazad. Semiconvergence and relaxation parameters for projected SIRT algorithms. SIAM J. Sci. Comput., (2012)34: A2000-A2017. [17] T. Elfving, T. Nikazad, and P. C. Hansen. Semi-convergence properties of Kaczmarz’s method. Inverse Problems, (2014)30: 055007. [18] Y. C. Eldar and D. Needell. Acceleration of randomized Kaczmarz method via the Johnson-Lindenstrauss Lemma. Numer. Alg., (2011)58: 163-177. [19] R. Escalante and M. Raydan. Alternating Projection Methods. SIAM Philadelphia, 2011. [20] G. H. Golub and C. V. Loan. Matrix Computations, 4th ed., Johns Hopkins University Press, Baltimore, MD, 2013. arik. Randomized iterative methods for linear systems. SIAM J. Matrix Anal. Appl., (2015)36: [21] R. M. Gower and P. Richt´ 1660-1690. [22] P. C. Hansen and M. S. Hansen. AIR Tools–A MATLAB package of algebraic iterative reconstruction methods. J. Comput. Appl. Math., (2012)236: 2167-2178. [23] P. C. Hansen and J. S. Jørgensen. AIR Tools II: algebraic iterative reconstruction methods, improved implementation. Numer. Alg. , (2018)79:107-137. [24] G. T. Herman. Fundamentals of Computerized Tomography: Image Reconstruction from Projections. 2nd edn (New York: Springer), 2009. [25] Y. Jiao, B. Jin, and X. Lu. Preasymptotic convergence of randomized Kaczmarz method. Inverse Problems, (2017)33: 125012. aherte Aufl¨ osung von Systemen linearer Gleichungen, in Bulletin International de l’Acad´ emie Polon[26] S. Kaczmarz. Angen¨ aise des Sciences et des Lettres. Classes des Sciences Mathmatiques et Naturelles. Series A: Sciences Mathmatiques, 1937: 355-357. [27] L. Landweber. An iterative formula for Fredholm integral equations of the first kind. Amer. J. Math., (1951)73: 615-624. [28] D. Leventhal and A. S. Lewis. Randomized methods for linear constraints: convergence rates and conditioning. Math.
Journal Pre-proof 12
Projected randomized Kaczmarz method
Jo
urn
al
Pr e-
p ro
of
Oper. Res., (2010)35: 641-654. [29] J. Lin and S. Li. Convergence of projected Landweber iteration for matrix rank minimization. Appl. Comput. Harm. Anal., (2014)36: 316-325. [30] R. Moazenzadeh, et al., Coupling a firefly algorithm with support vector regression to predict evaporation in northern Iran. Engineering Applications of Computational Fluid Mechanics, (2018)12: 584-597. [31] D. Needell. Randomized Kaczmarz solver for noisy linear systems. BIT Numer. Math., (2010)50: 395-403. [32] T. Nikazad, M. Abbasi, and T. Elfving. Error minimizing relaxation strategies in Landweber and Kaczmarz type iterations. J. Inverse Ill-Posed Probl., (2017)25: 35-56. [33] C. Popa. Least-squares solution of overdetermined inconsistent linear systems using Kaczmarz’s relaxation. Int. J. Comput. Math., (1995)55: 79-89. [34] C. Popa. Characterization of the solutions set of inconsistent least-squares problems by an extended Kaczmarz algorithm. Korean J. Comput. Appl. Math.,(1999)6: 51-64. [35] Y. Saad. Iterative Methods for Sparse Linear Systems. 2nd edition, SIAM, Philadelphia, PA, USA, 2000. [36] S. Samadianfard, et al., Daily global solar radiation modeling using data driven techniques and empirical equations in a semi-arid climate. Engineering Applications of Computational Fluid Mechanics, 2019 (13): 142-157. [37] R. W. Schafer, R. M. Mersereau, and M. A. Richards. Constrained iterative restoration algorithms. Proc. IEEE, (1981)69: 432-450. [38] T. Strohmer and R. Vershynin. A randomized Kaczmarz algorithm for linear systems with exponential convergence. J. Fourier Anal. Appl., (2009)15: 262-278. [39] C. Vogel. Computational Methods for Inverse Problems. SIAM, Philadelphia, 2002. [40] C.L. Wu, et al.,Rainfall-runoff modeling using artificial neural network coupled with singular spectrum analysis. J. Hydrol., ( 2011)399: 394-409. [41] H. Xiang and L. Zhang. Randomized iterative methods with alternating projections. arXiv preprint arXiv: 1708.09845, 2017. [42] Z.M. Yaseen, et al., An enhanced extreme learning machine model for river flow forecasting: state-of-the-art, practical applications in water resource engineering area and future research direction. J. Hydrol., (2019)569: 387-408. [43] H. Zhang and L. Z. Cheng. Projected Landweber iteration for matrix completion. J. Comput. Appl. Math., (2010)235: 593-601. [44] J. J. Zhang. A new greedy Kaczmarz algorithm for the solution of very large linear systems. Appl. Math. Letters, (2019)91: 207-212. [45] A. Zouzias and N. M. Freris. Randomized extended Kaczmarz for solving least-squares. SIAM J. Matrix Anal. Appl., (2012)34: 773-793. [46] X. Zheng, J. Yang, L. Li, et al. Modified projected landweber super-resolution algorithms for passive millimeter wave imaging. International Conference on Communications, 2008.