# On Randomized Sampling Kaczmarz Method with Application in Compressed Sensing

Authors: Mei-Lan Sun; Chuan-Qing Gu; Peng-Fei Tang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

We propose a randomized sampling Kaczmarz algorithm for the solution of very large systems of linear equations by introducing a maximal sampling probability control criterion, which is aimed at grasping the largest entry of the absolute sampling residual vector at each iteration. This new method differs from the greedy randomized Kaczmarz algorithm, which needs not to compute the residual vector of the whole linear system to determine the working rows. Numerical experiments show that the proposed algorithm has the most significant effect when the selected row number, i.e, the size of samples, is equal to the logarithm of all rows. Finally, we extend the randomized sampling Kaczmarz to signal reconstruction problems in compressed sensing. Signal experiments show that the new extended algorithm is more effective than the randomized sparse Kaczmarz method for online compressed sensing.


---

## Body

## 1. Introduction

In this paper, we mainly consider the iterative solution of large-scale consistent linear systems of the form(1)Ax=b,where A∈ℂm×n, b∈ℂm, and x∈ℂn denotes the n-dimensional unknown vector. The Kaczmarz algorithm [1] is a typical row-action method [2–5] for solving such a linear system (1) and is known as algebraic reconstruction technique in the field of computerized tomography [6, 7] and image reconstruction [5, 8–10]; see also [11–13] for additional references. More specifically, the (t + 1)th iterative vector xt + 1 is generated by the following Kaczmarz updates:(2)xt+1=xt+bit−AitxtAit22Ait∗,t=0,1,2,…,where Ait denotes the itth row of the matrix A, (⋅)∗ denotes the conjugate transpose of the corresponding vector or matrix, bit denotes the itth entry of the vector b, and it = (t mod m) + 1. Thus, the Kaczmarz algorithm is very suitable for the solution of very large linear equations since the main computation in the procedure of this algorithm is an inner product; however, this algorithm sometimes converges very slowly, see [14, 15] and the references therein. In order to improve the convergence of this algorithm, in 2009, Strohmer and Vershynin [16] proposed the randomized Kaczmarz (RK) algorithm with an expected exponential rate of convergence by using the rows of the coefficient matrix A in a random manner rather than in a given order, i.e., the RK chooses row it∈1,2,…,m with probability(3)Prrow=it=Ait22AF2.Recently, Bai and Wu proposed a different but more effective probability criterion and, based on it, they constructed the greedy randomized Kaczmarz (GRK) algorithm for solving the linear system (1) in [17]. This probability criterion makes the corresponding GRK algorithm to converge significantly faster than the RK method in both theory and experiments. We refer to [15, 18, 19] for more details on convergence theory and algorithmic generalizations for the GRK algorithm.We observe that the most expensive component at thetth iterate of the GRK algorithm is computing the square norm of the residual vector (i.e., b−Axt22) of the linear system, and it is the main factor that affects the computing time. Consequently, in this paper, we first randomly select k indices of rows, say, i1,i2,…,ik⊆1,2,…,, of the coefficient matrix of the linear system (1) and then compute the largest absolute relative residual, i.e, maxbil−Ailxt/AF, il∈i1,i2,…,ik and construct a randomized sampling Kaczmarz (RSK) algorithm. Therefore, our method needs not to calculate the residuals of all rows as the GRK method does and can be expected to be more effective than GRK in terms of the computing time. In theory, we provide a proof that shows linear convergence in expectation to the unique least-norm solution of the consistent linear system (1) for our RSK algorithm. And in computations, we show that the RSK algorithm significantly outperforms the RK algorithm in terms of both iteration counts and computing times and costs much less CPU time than the GRK algorithm but slightly requires more iteration steps than GRK algorithm. In addition, we extend the RSK algorithm to signal reconstruction problems in compressed sensing [20, 21] and establish the corresponding randomized sparse sampling Kaczmarz (RaSSK) algorithm and show that RaSSK performs much better than the randomized sparse Kaczmarz (RaSK) algorithm, which is a special case of the Bregman iterative algorithm [22] for online compressed sensing.The remainder of the paper is organized as follows. In Section2, we review the RK and GRK algorithms. In Section 3, we give the RSK algorithm and establish its convergence theory. Numerical experiments and comparisons are shown in Section 4. In Section 5, we extend the RSK to signal reconstruction in compressed sensing and give some signal experiments in this section. Finally, we end this paper with a few remarks and conclusions in Section 6.


## 2. The RK and GRK Algorithms

In this section, we briefly introduce the RK and GRK algorithms for solving systems of linear equations. First, let us give some necessary explanations. In this paper, we uset to denote the expected value conditional on the first t iterations, that is,(4)t.=deft.i1,i2,…,it,where il (l = 1, 2, …, t) is the ilth row chosen at the lth iterate. Then, it is easy to get that t.=. (see, e.g., the work of Bai and Wu [17]).When the linear system (1) is consistent and its coefficient matrix A∈ℂm×nm≥n is of the full-column rank, Strohmer and Vershynin [16] proposed the randomized Kaczmarz algorithm, as given in Algorithm 1.Algorithm 1:The randomized Kaczmarz (RK) algorithm.
Require:
A, b, l, and x0.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Selectit∈0,1,…,m with probability Prrow=it=Ait22/AF2
(3)
Setxt+1=xt+bit−Aitxt/Ait22Ait∗
(4)
end forThe RK method is convergent to the unique least-norm solution of the consistent linear system (1) when the coefficient matrix A is of the full-column rank with m ≥ n [16], and later in [23], Ma, Needell, and Ramdas gave the same convergence rate of the RK algorithm, but for the case that m < n. Compared with the Kaczmarz method, we can find, according to probability criterion (3), that the order in which the RK method selects the action row is based on the probability corresponding to the size of each row norm. It can improve the convergence of the Kaczmarz method discussed in [1].In [17], Bai and Wu proposed a new probability criterion (i.e., Steps 2–5 of Algorithm 2) for the RK algorithm, and based on it, they constructed a greedy randomized Kaczmarz algorithm, as given in Algorithm 2, which converges faster than the RK algorithm.Algorithm 2:The greedy randomized Kaczmarz (GRK) algorithm.
Require:
A, b, l, and x0.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Compute
ε
t
=
1
/
2
1
/
b
−
A
x
t
2
2
max
1
≤
i
t
≤
m
b
i
t
−
A
i
t
x
t
2
/
A
i
t
2
2
+
1
/
A
F
2
(3)
Determine the index set of positive integers
U
t
=
i
t
b
i
t
−
A
i
t
x
t
2
≥
ε
t
b
−
A
x
t
2
2
A
i
t
2
2
(4)
Compute theith entry r˜ti of the vector r˜t according to
r
˜
t
i
=
b
i
−
A
i
x
t
,
if
i
∈
U
t
,
0
,
otherwise
(5)
Selectit ∈ Ut with probability Prrow=it=r˜ti2/r˜t22
(6)
Setxt+1=xt+bit−Aitxt/Ait22Ait∗
(7)
end forAccording to Algorithm2, we can find that, at the tth iterate, the corresponding residual vector is rt = b − Axt and, if rti>rtj,i,j∈1,2,…,m, then the GRK algorithm is to make the ith row be selected with a larger probability than the jth row. Actually, only some entries of rt whose relative residual is greater than a threshold are nonzero. Convergence property of the GRK method has been studied by Bai and Wu in [17]. In addition, from [17], we can find that if the GRK method is convergent, then it must converge to the least-norm solution x∗ = A†b of the linear system (1), where A† indicates the Moore–Penrose pseudoinverse of the matrix A.


## 3. A Randomized Sampling Kaczmarz Method

### 3.1. The Randomized Sampling Kaczmarz Method

In this part, a new algorithm will be given. Through the contents of the previous sections, we note that at each iteration of GRK, it will cost a lot of money to calculate the residual vectorrt. If the dimension of the matrix is too large relative to the memory of the computer, the performance of the GRK method will be greatly reduced. In response to this situation, we propose a new method to solve the linear system (1). The purpose of our method is to make full use of the residual information as the GRK method does, while avoiding calculating the residual rt at each iteration. Let A=aij∈ℂm×n in the linear system (1) with 1 ≤ i ≤ m. Our method is carried out as follows: first, according to uniform distribution, randomly select k rows i1, i2, …, ik from the m rows of the coefficient matrix A and put them in a set Vt=i1,i2,…,ik; second, calculate the relative residuals of all corresponding rows in the set Vt; then, select the row with the largest absolute relative residual value (i.e., maxbil−Ailxt/AF, il∈i1,i2,…,ik) to implement the Kaczmarz update in (2). The RSK algorithm proposed in this paper is expressed in Algorithm 3.Algorithm 3:The randomized sampling Kaczmarz (RSK) algorithm.
Require:
A, b, l, k and x0.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Randomly generate by uniform distributionk different elements from 1 to m, indicated by Vt=i1,i2,…,ik
(3)
Compute thelth entry rtl of the vector rt according to
r
t
l
=
b
i
l
−
A
i
l
x
t
/
A
F
,
if
i
l
∈
V
t
,
l
=
1,2
,
…
,
k
(4)
Forl = 1, 2, …, k, select the row according to maxrtl for il ∈ Vt and assume that the index il = j
(5)
Setxt+1=xt+bi−Ajxt/Aj22Aj∗
(6)
end forAccording to Algorithm3, we can see that the selection of the appropriate number of rows k plays an important role in the RSK algorithm. More specifically, if the number of rows k is too small, it is difficult to ensure the efficiency of selecting iterative action rows in each iteration, and, if the number of rows k is too large, then the RSK method will have higher requirement for the memory of the computer, especially when the matrix dimension is large enough. In the following numerical examples, we will discuss the selection of parameter k in the RSK method.In fact, according to Algorithms1–3, the RSK method will cost (2k + 2) n + 1 + k flopping operations (flops) at each iteration, while the GRK method and the RK method will cost 7m + 2 (n + 1) and 4n + 1 flops, respectively [17]. We list the flops of these three algorithms at each iteration in Table 1. Obviously, when the row number m is extremely large, the RSK method costs less flops than the GRK method. When k = 1, the RSK method degenerates into the uniform sampling RK method.Table 1
The computational complexity of RK, GRK, and RSK algorithms at each iteration.
Algorithm
RK
GRK
RSK
Flopping operations (flops)
4n + 1
7m + 2(n + 1)
(2k + 2)n + 1 + k


### 3.2. Convergence Analysis of the RSK Algorithm

In this section, we pay attention to the convergence property of our algorithm. The following theorem guarantees the convergence of the RSK method.Theorem 1.
Letx∗ = A†b be the solution of the consistent linear system (1) and x0 ∈ R(A∗), where R(A∗) denotes the column space of the matrix A∗. Then, the iteration sequence xtt=0∞, generated by the RSK algorithm, converges to the solution x∗ in expectation. Moreover, the solution error in expectation for the iteration sequence xtt=0∞ obeys(5)xt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22,xt−x∗22≤1−λminA∗Am2×maxjAj22tx0−x∗22,t=0,1,2,…,where λminA∗A is the smallest positive eigenvalue of A∗A.Proof.
Update rule of the RSK algorithm yields(6)xt+1−xt=bit−AitxtAit22Ait∗,which indicates that xt+1 − xt is parallel to Ait∗. Using Aitx∗=bit, we compute(7)Aitxt+1−x∗=Aitxt−x∗+bit−AitxtAit22Ait∗=Aitxt−x∗+bit−Aitxt=bit−Aitx∗=0.
In other words,xt+1 − x∗ is orthogonal to Ait. Consequently, by the Pythagorean theorem, we obtain(8)xt+1−x∗22=xt−x∗22−xt+1−xt22.
Based on this equality, we have(9)txt+1−x∗22=xt−x∗22−txt+1−xt22=xt−x∗22−tbj−AjxtAj22Aj∗22=xt−x∗22−tbj−Ajxt|2Aj22.
Note that the population mean can be expressed as the expectation of the sample mean, thus,(10)tbj−Ajxt2Aj2=1ktkbj−Ajxt2Aj2≥1kt∑i=i1itbi−Aixt2Ai2=1m∑i=1mbi−Aixt2Ai2.
So, we upper bound(11)txt+1−x∗22≤xt−x∗22−1m2∑i=1mbi−Aixt|2maxjAj22=xt−x∗22−1m2Axt−x∗22maxjAj22=xt−x∗221−1m2Axt−x∗22maxjAj22xt−x∗22.
By the assumption thatx0 ∈ R(A∗), we have from the RSK algorithm that xt also does for each t. Since x∗ = A†b ∈ R(A∗), we obtain xt − x∗ ∈ R(A∗), which implies that(12)Axt−x∗22≥λminA∗Axt−x∗22.Hence, it holds that(13)txt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22.
In addition, by taking full expectation on both sides of (13) and using the law of iterated expectation k⋅=⋅, we immediately get(14)xt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22.
It follows from induction on the iteration indext, and we finally obtain that(15)xt−x∗22<1−λminA∗Am2×maxjAj22tx0−x∗22.Remark 1.
Actually, compared with the GRK algorithm in [17], the RSK algorithm cannot improve the convergence for the number of iterations for the linear system (1). However, when the row number m of the coefficient matrix A is much larger than k, in each iteration of the GRK algorithm, we must calculate the residuals of all m rows, while in the RSK algorithm, we just calculate the residual of k rows. This explains how and why the RSK algorithm uses less calculation time to achieve convergence accuracy when compared with the RK and GRK algorithms.


#### 3.2.1. Sensitivity Analysis

The essence of the RaSK algorithm is to solve the least square problem. Consider least square problem (16), given matrix A ∈ Rm×n, observation vector b ∈ Rm, and x ∈ Rn is the vector to be solved and satisfies(16)minAx−b2.LetA be accurate and vector b be the measured data, and then we consider the sensitivity of the solution x of problem (16) with respect to each component of data bj. Let the sensitivity of x with respect to bj be(17)Sij=∂xi∂bj.According to articles [24, 25], we have the following conclusions.Theorem 2.
In the least square problem (16), the sensitivity of x with respect to bj is(18)Sij=∂xi∂bj=∑k=1nuikλk∑l=1i≠knvlkbl+vjk,where uik and vlk are the elements of orthogonal matrix U and V, respectively, and it satisfies singular value decomposition formula A = VΛUT with Λ = diag (λ1, λ2, …) being the singular value of matrix A.Proof.
Without losing generality, let rank (A) = n(m ≥ n), and according to the singular value decomposition formula of the matrix, we have(19)A=VΛUT,where the orthogonal matrix V ∈ Rm×m, U ∈ Rn×n and Λ ∈ Rm×n.
LetΛ¯=diagλ1,λ2,…,λn and λ1 ≥ λ2 ≥…≥ λn, then we have(20)Λ=Λ¯0.
For anyx and b, let(21)α=UTx,β=VTb,then we can get(22)Ax−b=VλUTUα−Vβ=VΛ¯0α−β.
Letβ=β¯1/β¯2, where β¯1∈Rnandβ¯2∈Rm−n. We have(23)Ax−b=VΛ¯α−β¯1−β¯2,Ax−b22=Λ¯α−β¯1−β¯222=Λ¯α−β¯122+β¯222.
Sox should satisfy(24)Λ¯α=β¯1,Λ¯UTx=β¯1.
Thus, we can get(25)x=Λ¯UT−1β¯1=U−TΛ¯−1β¯1=UΛ¯−1β¯1,xi=∑k=1n1λkuikβk.
According toβ = VTb, we have(26)βk=∑l=1nvlkbl.
Thus, we can get(27)xi=∑k=1n1λkuik∑l=1nvlkbl=∑k=1n∑l=1nuikvkblλk=∑l=1n∑k=1nuikvlkλkbl,
So we have(28)∂xi∂bj=∑k=1nuikλk∑l=1l≠knvlkbl+vjk.
In fact,Sij reflects the rate of change of x with respect to data bj, which depends on the singular value of matrix A and its decomposition orthogonal matrix U, V. In particular, when m = n, the least square problem is reduced to solve the system of equations.


## 3.1. The Randomized Sampling Kaczmarz Method

In this part, a new algorithm will be given. Through the contents of the previous sections, we note that at each iteration of GRK, it will cost a lot of money to calculate the residual vectorrt. If the dimension of the matrix is too large relative to the memory of the computer, the performance of the GRK method will be greatly reduced. In response to this situation, we propose a new method to solve the linear system (1). The purpose of our method is to make full use of the residual information as the GRK method does, while avoiding calculating the residual rt at each iteration. Let A=aij∈ℂm×n in the linear system (1) with 1 ≤ i ≤ m. Our method is carried out as follows: first, according to uniform distribution, randomly select k rows i1, i2, …, ik from the m rows of the coefficient matrix A and put them in a set Vt=i1,i2,…,ik; second, calculate the relative residuals of all corresponding rows in the set Vt; then, select the row with the largest absolute relative residual value (i.e., maxbil−Ailxt/AF, il∈i1,i2,…,ik) to implement the Kaczmarz update in (2). The RSK algorithm proposed in this paper is expressed in Algorithm 3.Algorithm 3:The randomized sampling Kaczmarz (RSK) algorithm.
Require:
A, b, l, k and x0.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Randomly generate by uniform distributionk different elements from 1 to m, indicated by Vt=i1,i2,…,ik
(3)
Compute thelth entry rtl of the vector rt according to
r
t
l
=
b
i
l
−
A
i
l
x
t
/
A
F
,
if
i
l
∈
V
t
,
l
=
1,2
,
…
,
k
(4)
Forl = 1, 2, …, k, select the row according to maxrtl for il ∈ Vt and assume that the index il = j
(5)
Setxt+1=xt+bi−Ajxt/Aj22Aj∗
(6)
end forAccording to Algorithm3, we can see that the selection of the appropriate number of rows k plays an important role in the RSK algorithm. More specifically, if the number of rows k is too small, it is difficult to ensure the efficiency of selecting iterative action rows in each iteration, and, if the number of rows k is too large, then the RSK method will have higher requirement for the memory of the computer, especially when the matrix dimension is large enough. In the following numerical examples, we will discuss the selection of parameter k in the RSK method.In fact, according to Algorithms1–3, the RSK method will cost (2k + 2) n + 1 + k flopping operations (flops) at each iteration, while the GRK method and the RK method will cost 7m + 2 (n + 1) and 4n + 1 flops, respectively [17]. We list the flops of these three algorithms at each iteration in Table 1. Obviously, when the row number m is extremely large, the RSK method costs less flops than the GRK method. When k = 1, the RSK method degenerates into the uniform sampling RK method.Table 1
The computational complexity of RK, GRK, and RSK algorithms at each iteration.
Algorithm
RK
GRK
RSK
Flopping operations (flops)
4n + 1
7m + 2(n + 1)
(2k + 2)n + 1 + k


## 3.2. Convergence Analysis of the RSK Algorithm

In this section, we pay attention to the convergence property of our algorithm. The following theorem guarantees the convergence of the RSK method.Theorem 1.
Letx∗ = A†b be the solution of the consistent linear system (1) and x0 ∈ R(A∗), where R(A∗) denotes the column space of the matrix A∗. Then, the iteration sequence xtt=0∞, generated by the RSK algorithm, converges to the solution x∗ in expectation. Moreover, the solution error in expectation for the iteration sequence xtt=0∞ obeys(5)xt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22,xt−x∗22≤1−λminA∗Am2×maxjAj22tx0−x∗22,t=0,1,2,…,where λminA∗A is the smallest positive eigenvalue of A∗A.Proof.
Update rule of the RSK algorithm yields(6)xt+1−xt=bit−AitxtAit22Ait∗,which indicates that xt+1 − xt is parallel to Ait∗. Using Aitx∗=bit, we compute(7)Aitxt+1−x∗=Aitxt−x∗+bit−AitxtAit22Ait∗=Aitxt−x∗+bit−Aitxt=bit−Aitx∗=0.
In other words,xt+1 − x∗ is orthogonal to Ait. Consequently, by the Pythagorean theorem, we obtain(8)xt+1−x∗22=xt−x∗22−xt+1−xt22.
Based on this equality, we have(9)txt+1−x∗22=xt−x∗22−txt+1−xt22=xt−x∗22−tbj−AjxtAj22Aj∗22=xt−x∗22−tbj−Ajxt|2Aj22.
Note that the population mean can be expressed as the expectation of the sample mean, thus,(10)tbj−Ajxt2Aj2=1ktkbj−Ajxt2Aj2≥1kt∑i=i1itbi−Aixt2Ai2=1m∑i=1mbi−Aixt2Ai2.
So, we upper bound(11)txt+1−x∗22≤xt−x∗22−1m2∑i=1mbi−Aixt|2maxjAj22=xt−x∗22−1m2Axt−x∗22maxjAj22=xt−x∗221−1m2Axt−x∗22maxjAj22xt−x∗22.
By the assumption thatx0 ∈ R(A∗), we have from the RSK algorithm that xt also does for each t. Since x∗ = A†b ∈ R(A∗), we obtain xt − x∗ ∈ R(A∗), which implies that(12)Axt−x∗22≥λminA∗Axt−x∗22.Hence, it holds that(13)txt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22.
In addition, by taking full expectation on both sides of (13) and using the law of iterated expectation k⋅=⋅, we immediately get(14)xt+1−x∗22≤1−λminA∗Am2×maxjAj22xt−x∗22.
It follows from induction on the iteration indext, and we finally obtain that(15)xt−x∗22<1−λminA∗Am2×maxjAj22tx0−x∗22.Remark 1.
Actually, compared with the GRK algorithm in [17], the RSK algorithm cannot improve the convergence for the number of iterations for the linear system (1). However, when the row number m of the coefficient matrix A is much larger than k, in each iteration of the GRK algorithm, we must calculate the residuals of all m rows, while in the RSK algorithm, we just calculate the residual of k rows. This explains how and why the RSK algorithm uses less calculation time to achieve convergence accuracy when compared with the RK and GRK algorithms.


### 3.2.1. Sensitivity Analysis

The essence of the RaSK algorithm is to solve the least square problem. Consider least square problem (16), given matrix A ∈ Rm×n, observation vector b ∈ Rm, and x ∈ Rn is the vector to be solved and satisfies(16)minAx−b2.LetA be accurate and vector b be the measured data, and then we consider the sensitivity of the solution x of problem (16) with respect to each component of data bj. Let the sensitivity of x with respect to bj be(17)Sij=∂xi∂bj.According to articles [24, 25], we have the following conclusions.Theorem 2.
In the least square problem (16), the sensitivity of x with respect to bj is(18)Sij=∂xi∂bj=∑k=1nuikλk∑l=1i≠knvlkbl+vjk,where uik and vlk are the elements of orthogonal matrix U and V, respectively, and it satisfies singular value decomposition formula A = VΛUT with Λ = diag (λ1, λ2, …) being the singular value of matrix A.Proof.
Without losing generality, let rank (A) = n(m ≥ n), and according to the singular value decomposition formula of the matrix, we have(19)A=VΛUT,where the orthogonal matrix V ∈ Rm×m, U ∈ Rn×n and Λ ∈ Rm×n.
LetΛ¯=diagλ1,λ2,…,λn and λ1 ≥ λ2 ≥…≥ λn, then we have(20)Λ=Λ¯0.
For anyx and b, let(21)α=UTx,β=VTb,then we can get(22)Ax−b=VλUTUα−Vβ=VΛ¯0α−β.
Letβ=β¯1/β¯2, where β¯1∈Rnandβ¯2∈Rm−n. We have(23)Ax−b=VΛ¯α−β¯1−β¯2,Ax−b22=Λ¯α−β¯1−β¯222=Λ¯α−β¯122+β¯222.
Sox should satisfy(24)Λ¯α=β¯1,Λ¯UTx=β¯1.
Thus, we can get(25)x=Λ¯UT−1β¯1=U−TΛ¯−1β¯1=UΛ¯−1β¯1,xi=∑k=1n1λkuikβk.
According toβ = VTb, we have(26)βk=∑l=1nvlkbl.
Thus, we can get(27)xi=∑k=1n1λkuik∑l=1nvlkbl=∑k=1n∑l=1nuikvkblλk=∑l=1n∑k=1nuikvlkλkbl,
So we have(28)∂xi∂bj=∑k=1nuikλk∑l=1l≠knvlkbl+vjk.
In fact,Sij reflects the rate of change of x with respect to data bj, which depends on the singular value of matrix A and its decomposition orthogonal matrix U, V. In particular, when m = n, the least square problem is reduced to solve the system of equations.


## 3.2.1. Sensitivity Analysis

The essence of the RaSK algorithm is to solve the least square problem. Consider least square problem (16), given matrix A ∈ Rm×n, observation vector b ∈ Rm, and x ∈ Rn is the vector to be solved and satisfies(16)minAx−b2.LetA be accurate and vector b be the measured data, and then we consider the sensitivity of the solution x of problem (16) with respect to each component of data bj. Let the sensitivity of x with respect to bj be(17)Sij=∂xi∂bj.According to articles [24, 25], we have the following conclusions.Theorem 2.
In the least square problem (16), the sensitivity of x with respect to bj is(18)Sij=∂xi∂bj=∑k=1nuikλk∑l=1i≠knvlkbl+vjk,where uik and vlk are the elements of orthogonal matrix U and V, respectively, and it satisfies singular value decomposition formula A = VΛUT with Λ = diag (λ1, λ2, …) being the singular value of matrix A.Proof.
Without losing generality, let rank (A) = n(m ≥ n), and according to the singular value decomposition formula of the matrix, we have(19)A=VΛUT,where the orthogonal matrix V ∈ Rm×m, U ∈ Rn×n and Λ ∈ Rm×n.
LetΛ¯=diagλ1,λ2,…,λn and λ1 ≥ λ2 ≥…≥ λn, then we have(20)Λ=Λ¯0.
For anyx and b, let(21)α=UTx,β=VTb,then we can get(22)Ax−b=VλUTUα−Vβ=VΛ¯0α−β.
Letβ=β¯1/β¯2, where β¯1∈Rnandβ¯2∈Rm−n. We have(23)Ax−b=VΛ¯α−β¯1−β¯2,Ax−b22=Λ¯α−β¯1−β¯222=Λ¯α−β¯122+β¯222.
Sox should satisfy(24)Λ¯α=β¯1,Λ¯UTx=β¯1.
Thus, we can get(25)x=Λ¯UT−1β¯1=U−TΛ¯−1β¯1=UΛ¯−1β¯1,xi=∑k=1n1λkuikβk.
According toβ = VTb, we have(26)βk=∑l=1nvlkbl.
Thus, we can get(27)xi=∑k=1n1λkuik∑l=1nvlkbl=∑k=1n∑l=1nuikvkblλk=∑l=1n∑k=1nuikvlkλkbl,
So we have(28)∂xi∂bj=∑k=1nuikλk∑l=1l≠knvlkbl+vjk.
In fact,Sij reflects the rate of change of x with respect to data bj, which depends on the singular value of matrix A and its decomposition orthogonal matrix U, V. In particular, when m = n, the least square problem is reduced to solve the system of equations.


## 4. Numerical Experiments

### 4.1. Algorithm Comparison

In this section, we use RK, GRK, and RSK algorithms to solve equations (1) with different matrices A. The numerical behaviors of these algorithms are tested and evaluated in terms of the computing time in seconds (denoted as “CPU”) and the number of iteration steps (denoted as “IT”), and the CPU and IT mean the medians of the CPU time and iteration steps with respect to 50 times repeated runs of the corresponding method. In addition, we also report the speed-up of RSK (k) (k represents the number of rows selected from the coefficient matrix A in the RSK algorithm) against GRK, which is defined as(29)speed−up=CPUGRKCPURSKk.All algorithms started from the initial vectorx0 = 0 and terminated once the relative solution error (RSE), defined as(30)RSE=xt−x∗22x∗22,at the current iterate xt, satisfies RSE < 10−6, or the number of iteration steps exceeds 200,000. The latter is given a label “−” in the numerical tables. For part of the tested matrices, we give their Euclidean condition number, denoted by cond(A), and their density is defined as(31)density=number of nonzeros of anm×nmatrixm×n.All experiments are carried out using MATLAB (R2015b) on a personal laptop with 2.5 GHz (Intel(R) Core (TM) CPU i5-7300HQ), 8.00 GB memory, and Windows operating system (Windows 10).Example 1.
In this example, test matricesA are randomly generated by the MATLAB function randn, which produces independent entries subject to the standard normal distribution N (0,1). In our implementations, the random vectors x˜∈ℝn is randomly generated by using the MATLAB function randn, and b∈ℝm is taken to be Ax˜; in addition, the solution vectors x∗∈ℝn is taken to be pinv(A)b.
In Tables2 and 3, we report iteration counts and CPU times for RK, GRK, RSK (2), and RSK (7) when the linear system (1) is consistent. As the results in Table 2 show that the RSK (7) vastly outperform both the RK and GRK in terms of CPU times with significant speed-ups when the corresponding linear system is fat (i.e., m < n) and the minimum speed-up is 1.22 and the maximum reaches 2.08. From Table 3, we see that the CPU times of RSK (10) are considerably less than those of RK and GRK when the corresponding linear system is thin (i.e., m > n), with the speed-up being at least 2.32 and at most attaining 3.49. In addition, we can find that the “fatter” the matrix is, the RSK algorithm shows less advantages, and the “thinner” the matrix is, the RSK algorithm shows more advantages. It is in line with our intuition because if the row number m is extremely large, the RSK algorithm can reduce more computational complexity, for the RSK algorithm is independent of m while the GRK algorithm is not.Table 2
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with m = 100 and 200 and different n.
m × n
100 × 1000
100 × 3000
100 × 5000
200 × 1000
200 × 3000
200 × 5000
RK
IT
1127.3
906.0
825.6
3109.3
2137
2049.2
CPU
0.1335
0.1572
0.1907
0.3740
0.4142
0.5292
GRK
IT
319.4
232.1
212.5
932.6
552.6
481.4
CPU
0.0496
0.0626
0.0896
0.1574
0.2127
0.4336
RSK (2)
IT
739.9
595.4
548.0
2005.4
1374.4
1257.8
CPU
0.0396
0.0615
0.1043
0.1092
0.1907
0.3539
Speed-up
1.25
1.02
0.86
1.44
1.12
1.23
RSK (7)
IT
369.7
275.6
254.9
1060.8
649.2
583.8
CPU
0.0259
0.0456
0.0733
0.0757
0.1364
0.2564
Speed-up
1.92
1.37
1.22
2.08
1.56
1.69Table 3
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with n = 100 and 200 and different m.
m × n
1000 × 100
3000 × 100
5000 × 100
1000 × 200
3000 × 200
5000 × 200
RK
IT
1439.2
1329.5
1284.3
3550.4
2794.0
2655.7
CPU
0.1760
0.2407
0.3204
0.4416
0.5147
0.6745
GRK
IT
223.8
175.1
161.9
608.5
375.0
335.9
CPU
0.0377
0.0537
0.0809
0.1101
0.1553
0.2692
RSK (2)
IT
928.2
860.6
844.0
2246.0
1799.9
1733.7
CPU
0.0320
0.0365
0.0437
0.0803
0.0863
0.1057
Speed-up
1.18
1.47
1.85
1.37
1.80
2.55
RSK (10)
IT
399.6
372.3
368.5
978.1
773.8
753.5
CPU
0.0157
0.0231
0.0308
0.0460
0.0586
0.0771
Speed-up
2.40
2.32
2.63
2.39
2.65
3.49Example 2.
In this example, we select full-rank sparse matrices from [26], which originate in different applications such as linear programming, combinatorial optimization, DNA electrophoresis model, and Pajek or world city network. They possess certain structures and properties, such as square (m = n) (e.g., cage5), thin (m > n) (e.g., WorldCities), or fat (m < n) (e.g., bibd_16_8, and refine). There are also rank-deficient sparse matrices from [26], which come from applications like combinatorial problem and Pajek or world soccer network, such as square (e.g., football) and thin (e.g., relat6 and flower_5_1).
In Table4, we list the numbers of iteration steps and the computing times for RK, GRK, RSK (2), and RSK (5) algorithms. According to the test results in Table 4, we can see that the RSK algorithm performs better than the GRK algorithm in terms of CPU, even if the matrix is square or not, full-rank or not, and sparse or not and the condition number is infinite or not. More specifically, in terms of the RSK (2) algorithm, the speed-up is at least 1.59 and the biggest is 2.76; and in terms of the RSK (5) algorithm, the speed-up is at least 1.84 and the biggest even attains 4.45.Table 4
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with different m and n.
Name
refine
cage5
bibd−16−8
WorldCities
flower−5−1
relat6
football
m × n
29 × 62
37 × 37
120 × 12870
315 × 100
211 × 201
2340 × 157
35 × 35
Density
8.51%
17.02%
23.33%
23.87%
1.42%
2.21%
9.63%
Cond(A)
66.67
15.42
9.54
66.00
Inf
Inf
Inf
RK
IT
25963.3
10046.1
3058.5
24641
43611
9027.0
—
CPU
2.6477
1.0146
2.3381
2.5843
4.5635
1.4171
—
GRK
IT
689.1
556.1
1060.9
8925.8
9034.8
1672.9
102331
CPU
0.0858
0.0676
1.8835
1.3617
1.2300
0.4332
11.8878
RSK (2)
IT
1046.3
1173.5
1932.1
19249.1
17504.2
7578.5
189040.5
CPU
0.0311
0.0405
1.1868
0.5917
0.5764
0.2641
5.1743
Speed-up
2.76
1.67
1.59
2.30
2.13
1.64
2.30
RSK (5)
IT
632.4
614.3
1182.9
11435.0
11157.3
3808.3
127467.4
CPU
0.0193
0.0259
1.0222
0.3645
0.3762
0.1535
3.4984
Speed-up
4.45
2.61
1.84
3.74
3.27
2.82
3.40


### 4.2. The Choice of Parameterk

There is no doubt that the value of parameterk makes a huge difference in the performance of the RSK algorithm. In this section, we will have a tentative discussion on k and the matrix row number m. We simulate the relationship between the CPU of the RSK algorithm (with RSE < 10−6) and the size of k under different A. The simulation results are shown in Figure 1.Figure 1
Pictures of CPU for RSK versusk for a random matrix with different m and n. (a) Simulation results for m = 50, n = 1000, and log2 (m)≈5.6. (b) Simulation results for m = 1000, n = 50, and log2 (m)≈10. (c) Simulation results for m = 100, n = 1000, and log2 (m)≈6.6. (d) Simulation results for m = 2000, n = 0, and log2 (m)≈11. (e) Simulation results for m = 150, n = 1000, and log2 (m)≈7.2. (f) Simulation results for m = 3000, n = 50, and log2 (m)≈12.
(a)
(b)
(c)
(d)
(e)
(f)Through a large number of our numerical experiments, we find thatk = [log2(m)] may be a good choice.


## 4.1. Algorithm Comparison

In this section, we use RK, GRK, and RSK algorithms to solve equations (1) with different matrices A. The numerical behaviors of these algorithms are tested and evaluated in terms of the computing time in seconds (denoted as “CPU”) and the number of iteration steps (denoted as “IT”), and the CPU and IT mean the medians of the CPU time and iteration steps with respect to 50 times repeated runs of the corresponding method. In addition, we also report the speed-up of RSK (k) (k represents the number of rows selected from the coefficient matrix A in the RSK algorithm) against GRK, which is defined as(29)speed−up=CPUGRKCPURSKk.All algorithms started from the initial vectorx0 = 0 and terminated once the relative solution error (RSE), defined as(30)RSE=xt−x∗22x∗22,at the current iterate xt, satisfies RSE < 10−6, or the number of iteration steps exceeds 200,000. The latter is given a label “−” in the numerical tables. For part of the tested matrices, we give their Euclidean condition number, denoted by cond(A), and their density is defined as(31)density=number of nonzeros of anm×nmatrixm×n.All experiments are carried out using MATLAB (R2015b) on a personal laptop with 2.5 GHz (Intel(R) Core (TM) CPU i5-7300HQ), 8.00 GB memory, and Windows operating system (Windows 10).Example 1.
In this example, test matricesA are randomly generated by the MATLAB function randn, which produces independent entries subject to the standard normal distribution N (0,1). In our implementations, the random vectors x˜∈ℝn is randomly generated by using the MATLAB function randn, and b∈ℝm is taken to be Ax˜; in addition, the solution vectors x∗∈ℝn is taken to be pinv(A)b.
In Tables2 and 3, we report iteration counts and CPU times for RK, GRK, RSK (2), and RSK (7) when the linear system (1) is consistent. As the results in Table 2 show that the RSK (7) vastly outperform both the RK and GRK in terms of CPU times with significant speed-ups when the corresponding linear system is fat (i.e., m < n) and the minimum speed-up is 1.22 and the maximum reaches 2.08. From Table 3, we see that the CPU times of RSK (10) are considerably less than those of RK and GRK when the corresponding linear system is thin (i.e., m > n), with the speed-up being at least 2.32 and at most attaining 3.49. In addition, we can find that the “fatter” the matrix is, the RSK algorithm shows less advantages, and the “thinner” the matrix is, the RSK algorithm shows more advantages. It is in line with our intuition because if the row number m is extremely large, the RSK algorithm can reduce more computational complexity, for the RSK algorithm is independent of m while the GRK algorithm is not.Table 2
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with m = 100 and 200 and different n.
m × n
100 × 1000
100 × 3000
100 × 5000
200 × 1000
200 × 3000
200 × 5000
RK
IT
1127.3
906.0
825.6
3109.3
2137
2049.2
CPU
0.1335
0.1572
0.1907
0.3740
0.4142
0.5292
GRK
IT
319.4
232.1
212.5
932.6
552.6
481.4
CPU
0.0496
0.0626
0.0896
0.1574
0.2127
0.4336
RSK (2)
IT
739.9
595.4
548.0
2005.4
1374.4
1257.8
CPU
0.0396
0.0615
0.1043
0.1092
0.1907
0.3539
Speed-up
1.25
1.02
0.86
1.44
1.12
1.23
RSK (7)
IT
369.7
275.6
254.9
1060.8
649.2
583.8
CPU
0.0259
0.0456
0.0733
0.0757
0.1364
0.2564
Speed-up
1.92
1.37
1.22
2.08
1.56
1.69Table 3
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with n = 100 and 200 and different m.
m × n
1000 × 100
3000 × 100
5000 × 100
1000 × 200
3000 × 200
5000 × 200
RK
IT
1439.2
1329.5
1284.3
3550.4
2794.0
2655.7
CPU
0.1760
0.2407
0.3204
0.4416
0.5147
0.6745
GRK
IT
223.8
175.1
161.9
608.5
375.0
335.9
CPU
0.0377
0.0537
0.0809
0.1101
0.1553
0.2692
RSK (2)
IT
928.2
860.6
844.0
2246.0
1799.9
1733.7
CPU
0.0320
0.0365
0.0437
0.0803
0.0863
0.1057
Speed-up
1.18
1.47
1.85
1.37
1.80
2.55
RSK (10)
IT
399.6
372.3
368.5
978.1
773.8
753.5
CPU
0.0157
0.0231
0.0308
0.0460
0.0586
0.0771
Speed-up
2.40
2.32
2.63
2.39
2.65
3.49Example 2.
In this example, we select full-rank sparse matrices from [26], which originate in different applications such as linear programming, combinatorial optimization, DNA electrophoresis model, and Pajek or world city network. They possess certain structures and properties, such as square (m = n) (e.g., cage5), thin (m > n) (e.g., WorldCities), or fat (m < n) (e.g., bibd_16_8, and refine). There are also rank-deficient sparse matrices from [26], which come from applications like combinatorial problem and Pajek or world soccer network, such as square (e.g., football) and thin (e.g., relat6 and flower_5_1).
In Table4, we list the numbers of iteration steps and the computing times for RK, GRK, RSK (2), and RSK (5) algorithms. According to the test results in Table 4, we can see that the RSK algorithm performs better than the GRK algorithm in terms of CPU, even if the matrix is square or not, full-rank or not, and sparse or not and the condition number is infinite or not. More specifically, in terms of the RSK (2) algorithm, the speed-up is at least 1.59 and the biggest is 2.76; and in terms of the RSK (5) algorithm, the speed-up is at least 1.84 and the biggest even attains 4.45.Table 4
IT and CPU of RK, GRK, and RSK (k) for m × n matrices A with different m and n.
Name
refine
cage5
bibd−16−8
WorldCities
flower−5−1
relat6
football
m × n
29 × 62
37 × 37
120 × 12870
315 × 100
211 × 201
2340 × 157
35 × 35
Density
8.51%
17.02%
23.33%
23.87%
1.42%
2.21%
9.63%
Cond(A)
66.67
15.42
9.54
66.00
Inf
Inf
Inf
RK
IT
25963.3
10046.1
3058.5
24641
43611
9027.0
—
CPU
2.6477
1.0146
2.3381
2.5843
4.5635
1.4171
—
GRK
IT
689.1
556.1
1060.9
8925.8
9034.8
1672.9
102331
CPU
0.0858
0.0676
1.8835
1.3617
1.2300
0.4332
11.8878
RSK (2)
IT
1046.3
1173.5
1932.1
19249.1
17504.2
7578.5
189040.5
CPU
0.0311
0.0405
1.1868
0.5917
0.5764
0.2641
5.1743
Speed-up
2.76
1.67
1.59
2.30
2.13
1.64
2.30
RSK (5)
IT
632.4
614.3
1182.9
11435.0
11157.3
3808.3
127467.4
CPU
0.0193
0.0259
1.0222
0.3645
0.3762
0.1535
3.4984
Speed-up
4.45
2.61
1.84
3.74
3.27
2.82
3.40


## 4.2. The Choice of Parameterk

There is no doubt that the value of parameterk makes a huge difference in the performance of the RSK algorithm. In this section, we will have a tentative discussion on k and the matrix row number m. We simulate the relationship between the CPU of the RSK algorithm (with RSE < 10−6) and the size of k under different A. The simulation results are shown in Figure 1.Figure 1
Pictures of CPU for RSK versusk for a random matrix with different m and n. (a) Simulation results for m = 50, n = 1000, and log2 (m)≈5.6. (b) Simulation results for m = 1000, n = 50, and log2 (m)≈10. (c) Simulation results for m = 100, n = 1000, and log2 (m)≈6.6. (d) Simulation results for m = 2000, n = 0, and log2 (m)≈11. (e) Simulation results for m = 150, n = 1000, and log2 (m)≈7.2. (f) Simulation results for m = 3000, n = 50, and log2 (m)≈12.
(a)
(b)
(c)
(d)
(e)
(f)Through a large number of our numerical experiments, we find thatk = [log2(m)] may be a good choice.


## 5. Application of RSK Algorithm in Compressed Sensing

### 5.1. The RaSSK Algorithm

Consider the linear systemAx = b, where A∈ℝm×n,m<n, x is a sparse n dimension vector, and b is an m-dimension vector. By solving the following regularized basis pursuit problem(32)minx∈ℝnfx=λx1+12x22,s.t.Ax=b,we can find that the least Euclidean norm solution satisfies the sparsity requirement. In 2014, Lorenz et al. [22] proposed the RaSK algorithm for solving the regularized basis pursuit problem as given in Algorithm 4).Algorithm 4: The randomized sparse Kaczmarz (RaSK) algorithm [22].
Require:
x
0, A, b, and l.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Generateit randomly by Prrow=it=Ait22/AF2
(3)
Setxt+1/2=xt+bit−Aitxt/Ait22Ait∗
(4)
Setxt+1 = Sλ(xt+1/2)
(5)
end forShrink functionSλx=maxx|−λ,0⋅signx, where λ(>0) is set according to different applications and is to control the sparsity of the solution. Lorenz et al. proved that, for a consistent linear system Ax = b, the RaSK algorithm converges to the unique solution of the regularized basis pursuit problem in [22].Similar to the construction method of the RaSK algorithm, we can give the randomized sampling sparse Kaczmarz (RaSSK) algorithm which is listed in Algorithm5.Algorithm 5:The randomized sampling sparse Kaczmarz (RaSSK) algorithm.
Require:
x
0, A, b, l, k.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Randomly generate by uniform distributionk different elements from 1 to m, indicated by Vt=i1,i2,…,ik
(3)
Compute thelth entry rtl of the vector rt according to
r
t
l
=
b
i
l
−
A
i
l
x
t
/
A
F
,
if
i
l
∈
V
t
,
l
=
1,2
,
…
,
k
(4)
Forl = 1, 2, …, k, select the row according to maxrtl for il ∈ Vt and assume that the index il = j
(5)
Setxt+1/2=xt+bj−Ajxt/Aj22Aj∗
(6)
Setxt+1 = Sλ(xt+1/2)
(7)
end forThe proof of the convergence of RaSSK is similar to RaSK. It can be seen as a special case of the Bregman projections for split feasible problems (BPSFP) algorithm in [22]; if we change its feasibility question to “Find x∈∩k=1mAk” (i.e., the hyperplane formed by the rows of the matrix A) and define fx=λx1+1/2x22, then, the convergence can be easily obtained by Theorem 2.8 in [22].


### 5.2. Signal Experiments

In this section, we will show the efficiency of the RaSSK method by several numerical examples and compare it with the RaSK algorithm. We implement the numerical experiments, by MATLAB (R2015b) on a personal laptop with 2.5 GHz (Intel(R) Core(TM) CPU i5-7300HQ), 8.00 GB memory, and Windows operating system (Windows 10).Example 3.
In this example, the test signal is randomly generated with length 256 and limit its sparsity to 10, that is, only 10 nonzero coefficients. The signal is reconstructed by RaSK and RaSSK, respectively. The recovery signal map, relative error map, and relative residual map are given in Figure2. As shown in this figure, the parameter λ in the RaSSK method and the RaSK method is 50. We can see that both algorithms can reconstruct the original signal. It is worth mentioning that the RaSSK method is more efficient compared to the RaSK method.Figure 2
Recovery map, relative error map, and relative residual map in signal reconstruction by RaSK and RaSSK algorithms. (a) Signal recovery map by RaSK. (b) Signal recovery map by RaSSK. (c) Relative errorx−xreal2/xreal2. (d) Relative error Ax−b2/b2.
(a)
(b)
(c)
(d)Example 4.
In this example, we use the famous image lena with 64 × 64 and 128 × 128 pixels to test our algorithm. The MATLAB functionrandn generates an m × n random matrix which is independently distributed in the standard normal distribution N (0,1) as the observation matrix. For the sparse representation of the image, we use MATLAB function dwt, which is a wavelet transform process. The error is defined by(33)x−realx22realx22,and the stopping criterion of the algorithm is error < 0.1 or reaches the maximum number of iteration steps 50,000. The value of the parameter λ and the results are shown in Table 5. In this table, “CPU” denotes the computing time and “IT” denotes the number of iteration steps.
From Table5, it is easy to obtain that the RaSSK algorithm significantly performs better than the RaSK algorithm. It greatly reduces the number of iteration steps and the computing time. Meanwhile, from Figures 3 and 4, we can find that for human beings, there is almost no perception loss whenever some information is discarded.Table 5
Numerical results of RaSK and RaSSK for the 64 × 64 and 128 × 128 image to be reconstructed.
64 × 64
RaSK
RaSSK
128 × 128
RaSK
RaSSK
Observation matrix sizem × n
1000 × 2051
1000 × 2051
Observation matrix sizem × n
3000 × 8195
3000 × 8195
Original image sparsity
25.03%
25.03%
Original image sparsity
50.02%
50.02%
Sparse image sparsity
6.63%
6.63%
Sparse image sparsity
2.66%
2.66%
λ
50
50
λ
50
50
Iteration step
21483
47
Iteration step
50000
723
CPU
32.1677
0.0743
CPU
704.1457
10.6620Figure 3
The RaSK and RaSSK algorithms for a 64 × 64 reconstructed image. (a) is the original image, (b) is the image after sparse representation, (c) is the image obtained by the RaSK algorithm, and (d) is the image obtained by the RaSSK algorithm.
(a)
(b)
(c)
(d)Figure 4
The RaSK and RaSSK algorithms for a 128 × 128 reconstructed image. (a) is the original image, (b) is the image after sparse representation, (c) is the image obtained by the RaSK algorithm, and (d) is the image obtained by the RaSSK algorithm.
(a)
(b)
(c)
(d)


## 5.1. The RaSSK Algorithm

Consider the linear systemAx = b, where A∈ℝm×n,m<n, x is a sparse n dimension vector, and b is an m-dimension vector. By solving the following regularized basis pursuit problem(32)minx∈ℝnfx=λx1+12x22,s.t.Ax=b,we can find that the least Euclidean norm solution satisfies the sparsity requirement. In 2014, Lorenz et al. [22] proposed the RaSK algorithm for solving the regularized basis pursuit problem as given in Algorithm 4).Algorithm 4: The randomized sparse Kaczmarz (RaSK) algorithm [22].
Require:
x
0, A, b, and l.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Generateit randomly by Prrow=it=Ait22/AF2
(3)
Setxt+1/2=xt+bit−Aitxt/Ait22Ait∗
(4)
Setxt+1 = Sλ(xt+1/2)
(5)
end forShrink functionSλx=maxx|−λ,0⋅signx, where λ(>0) is set according to different applications and is to control the sparsity of the solution. Lorenz et al. proved that, for a consistent linear system Ax = b, the RaSK algorithm converges to the unique solution of the regularized basis pursuit problem in [22].Similar to the construction method of the RaSK algorithm, we can give the randomized sampling sparse Kaczmarz (RaSSK) algorithm which is listed in Algorithm5.Algorithm 5:The randomized sampling sparse Kaczmarz (RaSSK) algorithm.
Require:
x
0, A, b, l, k.
Ensure:
x
l.
(1)
for
t = 0, 1, …, ldo
(2)
Randomly generate by uniform distributionk different elements from 1 to m, indicated by Vt=i1,i2,…,ik
(3)
Compute thelth entry rtl of the vector rt according to
r
t
l
=
b
i
l
−
A
i
l
x
t
/
A
F
,
if
i
l
∈
V
t
,
l
=
1,2
,
…
,
k
(4)
Forl = 1, 2, …, k, select the row according to maxrtl for il ∈ Vt and assume that the index il = j
(5)
Setxt+1/2=xt+bj−Ajxt/Aj22Aj∗
(6)
Setxt+1 = Sλ(xt+1/2)
(7)
end forThe proof of the convergence of RaSSK is similar to RaSK. It can be seen as a special case of the Bregman projections for split feasible problems (BPSFP) algorithm in [22]; if we change its feasibility question to “Find x∈∩k=1mAk” (i.e., the hyperplane formed by the rows of the matrix A) and define fx=λx1+1/2x22, then, the convergence can be easily obtained by Theorem 2.8 in [22].


## 5.2. Signal Experiments

In this section, we will show the efficiency of the RaSSK method by several numerical examples and compare it with the RaSK algorithm. We implement the numerical experiments, by MATLAB (R2015b) on a personal laptop with 2.5 GHz (Intel(R) Core(TM) CPU i5-7300HQ), 8.00 GB memory, and Windows operating system (Windows 10).Example 3.
In this example, the test signal is randomly generated with length 256 and limit its sparsity to 10, that is, only 10 nonzero coefficients. The signal is reconstructed by RaSK and RaSSK, respectively. The recovery signal map, relative error map, and relative residual map are given in Figure2. As shown in this figure, the parameter λ in the RaSSK method and the RaSK method is 50. We can see that both algorithms can reconstruct the original signal. It is worth mentioning that the RaSSK method is more efficient compared to the RaSK method.Figure 2
Recovery map, relative error map, and relative residual map in signal reconstruction by RaSK and RaSSK algorithms. (a) Signal recovery map by RaSK. (b) Signal recovery map by RaSSK. (c) Relative errorx−xreal2/xreal2. (d) Relative error Ax−b2/b2.
(a)
(b)
(c)
(d)Example 4.
In this example, we use the famous image lena with 64 × 64 and 128 × 128 pixels to test our algorithm. The MATLAB functionrandn generates an m × n random matrix which is independently distributed in the standard normal distribution N (0,1) as the observation matrix. For the sparse representation of the image, we use MATLAB function dwt, which is a wavelet transform process. The error is defined by(33)x−realx22realx22,and the stopping criterion of the algorithm is error < 0.1 or reaches the maximum number of iteration steps 50,000. The value of the parameter λ and the results are shown in Table 5. In this table, “CPU” denotes the computing time and “IT” denotes the number of iteration steps.
From Table5, it is easy to obtain that the RaSSK algorithm significantly performs better than the RaSK algorithm. It greatly reduces the number of iteration steps and the computing time. Meanwhile, from Figures 3 and 4, we can find that for human beings, there is almost no perception loss whenever some information is discarded.Table 5
Numerical results of RaSK and RaSSK for the 64 × 64 and 128 × 128 image to be reconstructed.
64 × 64
RaSK
RaSSK
128 × 128
RaSK
RaSSK
Observation matrix sizem × n
1000 × 2051
1000 × 2051
Observation matrix sizem × n
3000 × 8195
3000 × 8195
Original image sparsity
25.03%
25.03%
Original image sparsity
50.02%
50.02%
Sparse image sparsity
6.63%
6.63%
Sparse image sparsity
2.66%
2.66%
λ
50
50
λ
50
50
Iteration step
21483
47
Iteration step
50000
723
CPU
32.1677
0.0743
CPU
704.1457
10.6620Figure 3
The RaSK and RaSSK algorithms for a 64 × 64 reconstructed image. (a) is the original image, (b) is the image after sparse representation, (c) is the image obtained by the RaSK algorithm, and (d) is the image obtained by the RaSSK algorithm.
(a)
(b)
(c)
(d)Figure 4
The RaSK and RaSSK algorithms for a 128 × 128 reconstructed image. (a) is the original image, (b) is the image after sparse representation, (c) is the image obtained by the RaSK algorithm, and (d) is the image obtained by the RaSSK algorithm.
(a)
(b)
(c)
(d)


## 6. Conclusion

Variants of the RK algorithm are effective iteration methods for large-scale linear systems. In this paper, based on the randomized Kaczmarz algorithm and the greedy randomized Kaczmarz algorithm, we propose a new algorithm which makes use of the residual information, while it need not calculate all the residuals. This algorithm converges faster than RK [16] and GRK [17] in experiments. Furthermore, after a large amount of numerical experiments, we recommend the parameter k to take [log2(m)]. As an application, we apply it to the signal reconstruction in compressed sensing. The experiments show that our algorithm has a good performance.



---

*Source: 7464212-2020-03-30.xml*