A two-step randomized Gauss-Seidel method for solving large-scale linear least squares problems

1.   Introduction

We consider the approximate solutions of a large linear least squares problem

where $A\in\mathbb{R}^{m\times n}$, $m > n$, is of full clolumn rank. $A$ and $b\in\mathbb{R}^m$ are known, $x\in\mathbb{R}^n$ is unknown to be determined. Under the condition that the linear system is consistent or inconsistent, people are interested in finding the unique least squares solution $x_* = A^{\dagger}b$, where $A^{\dagger}$ is the Moore-Penrose pseudoinverse of the matrix $A$ ($A^{\dagger} = (A^TA)^{-1}A^T$, where $A^T$ denotes the transpose of the matrix $A$). See also references ^{[1,2,3,4,5,6]}. As we know, the least squares problem arises widely in many fields such as tomography^[7,8,9], protein structure^[10], machine learning^[11], biological feature selection^[12], and so on^[13,14,15]. For solving (1.1), there are many direct methods, such as QR decomposition and singular value decomposition (SVD) ^[1,16]. However, for large-scale system matrices, these methods are too expensive because they consume a lot of memory space. So, some iterative methods are applied to solve large-scale linear least squares problems.

As one of the famous iterative algorithms, the Gauss-Seidel method ^[17] selects a coordinate $d_k$ and a step $\alpha_k = arg \min\limits_{\alpha\in\mathbb{R}}f(x_k+\alpha d_k)$ in each iteration. When $\alpha_k = A_{j_k}^T(b-Ax_k)/\left\|A_{j_k}\right\|^2_2$ and $d_k = e_{j_k}$ are accurately given, it generates the following iterative process:

where $j_k = (k\; mod\; n)+1$, $(\cdot)^T$ represents the transpose of a matrix or a vector, and $e_{j_k}$ represents the coordinate column vector, whose $j_k$-th entry is $1$ and zero otherwise. Inspired by the randomized Kaczmarz (RK) linear convergence characteristics of Strohmer and Vershynin^[18], Leventhal and Lewis^[19] obtained a randomized Gauss-Seidel (RGS) method, which is also known as randomized coordinate descent (RCD) method to solve (1.1). The RGS selects the update column index $j_k$ according to the appropriate probability. Theoretical analysis shows that RGS converges linearly to the unique least squares solution $x_* = A^{\dagger}b$. Many variants of RGS are receiving extensive attention recently due to its good performance. For example, the versions of block^[20,21], random greedy^[22,23,24]. Other versions, see literatures ^[4,17,25,26] and reference therein.

In 2018, Wu^[27] obtained a randomized block Gauss-Seidel (RBGS) method, which can significantly improve the convergence speed. However, the good column pavings for the RBGS is difficult to find when the column norms of the coefficient matrix fluctuate in a large range. In this paper, we propose a two-step randomimzed Gauss-Seidel method (TRGS), which does not need any columns pavings. The convergence of the TRGS algorithm is proved for (1.1) with the coefficient matrix of full rank.

This paper is organized as follows. In Section 2, some symbols and a lemma related to RGS are introduced. The convergence of RGS2 is analyzed. In Section 3, a convergent upper bound of TRGS is obtained theoretically. Several numerical experiments are reported to verify the feasibility of our proposed algorithms in Section 4. The conclusions of this paper are given in Section 5 and the proof of Theorem 2 is shown in appendix.

2.   The simplified two-step randomized Gauss-Seidel method

We first introduce the notations and definitions as follows. For the matrix $Q\in \mathbb{R}^{m\times n}$, $Q_{i}$ represents the $i$th column of the $Q$, $\lambda_{max}(Q^TQ)$ and $\lambda_{min}(Q^TQ)$ represent the maximum and minimum positive eigenvalues of the $Q^TQ$, respectively, $(\cdot)^T$ represents the transpose of a matrix or a vector, the $Q_{\tau_k}$ represents a submatrix of the $Q$ (where $\tau_k$ is a set of column indexes), $\|Q\|_2$ and $\|Q\|_F$ represent the spectral norm and Frobenius norm of the $Q$, respectively, and $\mathcal{R}(Q)$ is the image space of matrix $Q$. For a vector $p\in\mathbb{R}^n$ or $p\in\mathbb{R}^m$, $p_i$ is the ith component of $p$. For constant $c\in\mathbb{R}$, $[c]$ refers to the set consising of all positive integers not exceeding $c$. For the matrix in (1.1), assuming that every two columns are independent and identically distributed, the correlation coefficient parameters of the matrix $A$ are defined as follows

Then, we have $0\leq\delta\leq\Delta < 1$.

The randomized Gauss-Seidel method proposed by Leventhal and Lewis^[19] consists of two parts. The RGS determines a column index $j_k$ according to the probability $Pr(column = j_{k}) = \frac{\|A_{j_{k}}\|^2_2}{\|A\|^2_F}$ and then updates $x_{k+1}$ by (1.2).

The following results in ^[19] summarized an upper bound for the error of the solution in expectation on the convergence of RGS algorithm.

Lemma 1. Assume the least squares problems $(1.1)$ has tall and narrow coefficient matrix $A\in\mathbb{R}^{m\times n}$ with full rank. Let $x_* = A^{\dagger}b$ be a solution of $(1.1)$. Given an initial guess $x_0\in\mathbb{R}^n$, then the sequence $\{x_k\}_{0}^{\infty}$ generated by RGS algorithm converges linearly to the $x_*$ in expectation. Moreover, it satisfies

Now, we propose a simplified two-step randomized Gauss-Seidel method (RGS2) to solve (1.1). Algorithm 1 lists the RGS2 method, which mainly consists of two stages. The first stage is to select two different working columns by

where $j_{k_1}\in\{1, 2, \cdots, n\}$, $j_{k_2}\in\{1, 2, \cdots, n\}/\{j_{k_1}\}$. The second stage is to use (1.2) twice to update $x_k$.

If the coefficient matrix $A$ in (1.1) is tall and narrow with full column rank, the following result gives the convergence of the RGS2 algorithm.

Theorem 1. Assume the least squares problem $(1.1)$ has tall and full-rank coefficient matrix $A\in\mathbb{R}^{m\times n}$. Let $x_* = A^{\dagger}b$ be a solution of $(1.1)$. Given an initial guess $x_0\in\mathbb{R}^n$, then the iterative sequence $\{x_k\}_{0}^{\infty}$ generated by RGS2 algorithm converges linearly to the $x_*$ in expectation. Moreover, it satisfies

and

where $\tau_{max} = \|A\|^2_F-\min\limits_{p\in[n]}\|A_{p}\|^2_2$.

Proof. Let $\hat{y}_k$ and $\hat{x}_{k+1}$ be the first and the second iterative solutions of $x_k$ obtained by single-step continuous execution of Algorithm 1, respectively. The update process is divided into two steps as follows

Let $P_{j_{k_i}} = I_m-\frac{A_{j_{k_i}}A_{j_{k_i}}^T}{\|A_{j_{k_i}}\|^2_2}$ satisfy $P_{j_{k_i}}^2 = P_{j_{k_i}}, P_{j_{k_i}}^T = P_{j_{k_i}}$, where $i = 1, 2.$ Then, $P_{j_{k_i}}$ is the projection matrix, and

where the second equality follows from the normal equation $A^TAx_* = A^Tb$, whose $j_{k_2}$-th equation gives $A_{j_{k2}}^Tb = A_{j_{k_2}}^TAx_*$. Taking the expectation on the equality above, and by the expectation conditional upon $j_{k_1}$ (we fix the choice of $j_{k_1}$ and averagever the random index $j_{k_2}$), one can get

in which the last equation holds because

Note that

then $\|I_m- {AA^T\over \|A\|_F^2-\|A_{j_{k_1}}\|_2^2}\|_2\le 1-{\lambda_{\rm min}(A^TA)\over \tau_{\rm max}}$ and

Therefore,

(2.2) is obtained from the recurrence relation of (2.3) and the full expectation formula. This completed the proof.

3.   Two-step randomized Gauss-Seidel method

The minimum positive eigenvalue of $\lambda_{min}(A^TA)$ will become very small if the matrix $A$ has high correlation parameters^[28]. A weak bound of the convergence in Theorem 1 will appear. Under this condition, the angle of the unit coordinate direction as the search direction for two consecutive stages may be too small, which is the main reason for the slow convergence of RGS. Inspired by the work of Needell^[28] and Wu^[29], we iteratively update the solution by continuously seeking two more extensive directions, that is, determine the column pairs $(r, s)$ in advance. So a two-step randomized Gauss-Seidel algorithm is obtained. In a two-step randomized algorithm, which mainly consists of three steps as follows: First, we randomly select the two columns $r, s\in[n]$ according to the probability criterion, then estimate the optimal parameter $\lambda_{opt}$ for the first iteration with $y_k = x_k+\lambda_{opt}\frac{A_r^Tr_k}{\|A_r\|^2_2}e_r$, and finally perform the second iteration to update the iterative solution by $x_{k+1} = y_k+\frac{A_s^T(b-Ay_k)}{\|A_s\|^2_2}e_s$.

We consider the convergence of the TRGS algorithm. We first need the following result presented in ^[28].

Lemma 2. For any $\epsilon\in\mathbb{R}$, $\phi, \psi\in\mathbb{R}^m$, the minimizer of $\|\epsilon\phi+\psi\|^2_2$ is $\epsilon_{opt} = -\frac{ < \phi, \psi > }{\|\phi\|^2_2}$.

Now, we note that

then

where $\mu_k = \frac{(A_r)^T(A_s)}{\|A_r\|_2\|A_s\|_2}$. Obviously,

The selection of optimal parameter $\lambda_{opt}$ aims to minimize $\|Ax_{k+1}-b\|^2_2$. By Lemma 2, one can get the optimal value of $\lambda_k$, that is,

Substituting $\lambda_{opt}$ into (3.1), we obtain

So,

Set $r = j_{k_1}$ and $s = j_{k_2}$ with the optimal parameter $\lambda_{opt}$, this process can be described as

where

In order to simplify the computation, we rewrite the iterative process in Algorithm 2.

The following Lemmas 3 and 4 will be used to prove the convergence of TRGS algorithm.

Lemma 3. For the column vector $u_k$ defined above, $\|u_k\|^2_2 = 1-\mu_k^2$.

Proof. By the definition of $\mu_k$ in (3.1) and $u_k$ in (3.2), one can obtain that

Lemma 4. Define $\alpha_{s, t}$ and $\beta_{s, t}$ as

then, the existence of $\gamma\in\mathbb{R}$ subject to $(\vert \alpha_{s, t}\vert -\vert \beta_{s, t}\vert)^2\geq\gamma$.

Proof. It is known that $\delta = \min\limits_{s\neq t}\frac{\vert A_s^TA_t\vert }{\|A_s\|_2\|A_t\|_2}$ and $\Delta = \max\limits_{s\neq t}\frac{\vert A_s^TA_t\vert }{\|A_s\|_2\|A_t\|_2}$, Then,

When the coefficient matrix $A$ in (1.1) is tall and narrow with full column rank, the following result gives the convergence of the TRGS algorithm.

Theorem 2. Assume that the tall and narrow coefficient matrix $A\in R^{m\times n}$ has full column rank. Let $x_* = A^{\dagger}b$ be a solution of $(1.1)$. Given an initial guess $x_0\in\mathbb{R}^n$, then the sequence $\{x_k\}_{0}^{\infty}$ generated by TRGS algorithm converges linearly to the $x_*$ in expectation. Moreover, it satisfies

and

where $\tau_{max} = \max\limits_{p\in[n]}\{\|A\|^2_F-\|A_{p}\|^2_2\}$, $\tau_{min} = \min\limits_{q\in[n]}\{\|A\|^2_F-\|A_{q}\|^2_2\}$ and $\gamma = \min\begin{Bmatrix} \frac{\delta^2(1-\delta)}{1+\delta}, \frac{\Delta^2(1-\Delta)}{1+\Delta} \end{Bmatrix}$.

Proof. See Appendix 1.

Remark 1. We remark that the upper bound of the convergence of RGS method from Lemma 1 is

From Theorem 1, we remark the upper bound of the convergence of RGS2 method is

By Theorem 2, we can obtain an upper bound of the convergence of TRGS method

At each iteration, the TRGS method uses two columns of the matrix, while RGS utilizes only one. To be fair, we compare the upper bound on the convergence factor of one iteration of the TRGS method with that of two iterations, instead of one iteration, of the RGS method. Note that $0\leq\gamma < 1$ and $\tau_{max}\leq\|A\|^2_F$, we have $0\leq\gamma\tau_{min}/\tau_{max} < 1$ then $\Psi_{TRGS}\leq\Psi_{RGS2}\leq\Psi_{RGS}^2 < \Psi_{RGS}$. Especially, when the column correlation coefficient $\delta\text{ or }\Delta = 0$, then $\gamma = 0$. The RGS2 and TRGS methods have the same convergence factor in the upper bound.

Remark 2. If $\|A_{j}\|^2_2, \; j = 1, \cdots, n$, are precomputed, we discuss the computational cost of RGS, RGS2 and TRGS in each iteration step. The RGS costs $2m+2n+1$ flopping operations (flops), the RGS2 method needs $4m+4n+2$ flops, while the TRGS method requires $6m+4n+11$ flops.

4.   Numerical examples

In this section, we give several examples to show the efficiency of our TRGS method. We compare TRGS with RGS2 and RGS. In addition, randomized Kaczmarz (RK) in ^[18], randomized extended Kaczmarz (REK) in ^[2], partially randomized extended Kaczmarz (PREK) in ^[30], generalized two-subspace randomized Kaczmarz (GTRK) and two-subspace randomized extended Kaczmarz (TREK) in ^[29], as other iterative methods, are considered for solving consistent or inconsistent linear systems. All experiments are carried out with the Matlab 2020b on a computer with 3.00 GHz processing unit and 16 GB RAM.

We measure the efficiency of TRGS and other methods by the relative solution error

The initial vector is set as $x_0 = (0, 0, \cdots, 0)^T$ for all methods. When the set maximum number of iterations $k_{\text{max}} = 10^6 \text{ or RSE} < 10^{-6}$, we terminate the iteration process of each method. The '-' means that the number of iteration steps of the algorithm reaches $k_{\text{max}}$.

Example 4.1. In this example, for conherent matrix $A\in\mathbb{R}^{2000\times 100}$ in least-squares problem (1.1). The entries of the $A$ are the independent identically distributed uniform random variables in the interval $(t, 1)$, where the $t\in[0, 1]$. We remark the average column correlation index as follows

When $t$ increases from 0 to 1, the change of the $\bar{\mu}_k$ with $t$ and the relationship between $\Psi_{\text{RGS}}$, $\Psi_{\text{RGS}}^2$, $\Psi_{\text{RGS2}}$ and $\Psi_{\text{TRGS}}$ versus $t$ are plotted in Figure 1.

As can be seen from Figure 1, the left subfigure shows that as $t$ approaches 1, the average column correlation index $\bar{\mu}_k$ is highly correlated. From the right subfigure, the convergence upper bound of RGS2 and TRGS are always lower than RGS, and further we find that the theoretical upper bound of convergence of the TRGS will not exceed the RGS2.

Example 4.2 We use the RGS2 and TRGS methods to test the consistent least squares problem (1.1) with different size and compare them with the RK, RGS and the GTRK methods. The size of $A$ is $m\times n$ with $m = 10^3*k$ $(k = 1, 2, \cdots, 5)$ and $n = 50$. The entries of the matrix $A$ are the independent identically distributed uniform random variables in the interval $(t, 1)$. The vector $b = Ax_*$, where $x_*$ is generated randomly with the MATLAB function randn.

In Tables 1–3, we list the IT and the CPU(s) for the RK, RGS, RGS2, GTRK, and TRGS methods with $t = 0.1, 0.5\text{ and }0.8$, and the Euclidean condition number Cond(A) of the matrix is reported in each table. Figure 2 shows the plots of $m$ versus IT (left), and m versus CPU (right) of Algorithm 2 applied to solve (1.1) with different coefficient matrix $A$ listed in Tables 1–3, respectively.

From Tables 1–3, we can see that the TRGS method is better than other algorithms based on IT and CPU(s). We find that GTRK and TRGS are basically stable in both IT and CPU(s), while RK, RGS and RGS2 methods need more iterations and CPU time.

From Figure 2, it is not difficult to see that the curve of TRGS is much lower than that of RGS2 in terms of the IT and the CPU(s). In addition, RGS2 is sensitive to $t$, while TRGS is not affected by it. For example, in Figure 2, fix $m = 3000$, when $t = 0.1, \; 0.5 \text{ and } 0.8$, the IT of TRGS are basically steady at the level about $10^7$, while the IT of RGS2 are steady at the level $10^8, 10^{10} \text{ and } 10^{12}$, respectively. Similar results also appear in the CPU(s).

Example 4.3 In this example, we apply the RGS2 and TRGS methods to solve the inconsistent least squares problem (1.1) and compare them with the REK, RGS, PREK and TREK methods. The entries of the matrix $A$ are the independent identically distributed uniform random variables in the interval $(t, 1)$, and the vector $b = Ax_*+r$, where $x_*$ is one of the solutions of (1.1), which is generated randomly with the MATLAB function randn, and $r$ is a nonzero vector in the null space of $A^T$, which is generated by the MATLAB function null. The size of $A$ is $m\times n$ with $m = 10^3*k$ $(k = 1, 2, \cdots, 5)$, and $n = 100$.

Tables 4–6 list the iteration number (IT) and the CPU time when all methods stop and we also set $t = 0.1, 0.5, 0.8$ in each case. Figure 3 shows the plots of $m$ versus IT (left) and $m$ versus CPU (right) of TRGS and RGS2 applied to solve all linear systems (1.1) in Tables 4–6.

From Tables 4–6, we can see that the TRGS method is better than other algorithms based on IT and CPU(s). We also find that TREK and TRGS are basically stable in both IT and CPU(s), while the REK, RGS, PREK and RGS2 need more iterations and CPU time.

From Figure 3, the curves of TRGS and RGS2 show that RGS2 needs more iterations and CPU time to reach the stopping criterion. In addition, RGS2 is sensitive to $t$, while TRGS is not. For example, in Figure 3, fix $m = 3000$, when $t = 0.1, \; 0.5 \text{ and } 0.8$, the IT of TRGS are basically steady at the level about $10^7$, while the IT of RGS2 are steady at the level about $10^8, 10^{10} \text{ and } 10^{12}$ respectively. Similar results also appear in the CPU(s).

Example 4.4 In this example, we apply the TRGS method to solve the least squares problem (1.1) with the sparse coefficient matrix $A$ taken from the Florida sparse matrix collection in ^[31]. Especially, we select the tall and narrow sparse matrix $A$ with full column rank. Table 7 summarizes different sparse systems with density and condition number $Cond(A)$, where the density of a matrix $A$ means the ratio of the number of the nonzero elements of $A$ to the total number of the elements of $A$. Algorithm 2 is compared with the RK, RGS, GTRK, REK, PREK, TREK and RGS2 methods.

When the sparse least squares problem (1.1) is set to be consistent, the vector $b = Ax_*$, where $x_*$, one of the solutions of learst squares problem (1.1), is generated randomly with the MATLAB function randn. Table 8 lists the iteration number (IT) and the CPU time when RK, RGS, RGS2, GTRK and TRGS methods stop. Figure 4 shows the plot of RSE versus IT (left) and RSE versus CPU (right) of Algorithm 2 applied to solve (1.1) with the sparse coefficient matrix $A$ named "divorce".

When the sparse least squares problem (1.1) is set to be inconsistent, the $b = Ax_*+r$, where the $r$ is a nonzero vector in the null space of $A^T$. Due to the large dimension of $A$, the $r$ cannot be generated by the Matlab function null, but it can be generated by the projection vector $\check{r}$. The $\check{r}$ is generated by the MATLAB function randn and projected onto the null space of $A^T$. That is to say, $r = \check{r}-AA^{\dagger}\check{r}$, where $A^{\dagger}\check{r}$ is obtained by the Matlab function lsqminnorm. Table 9 lists the IT and the CPU(s) when REK, RGS, PREK, RGS2, TREK and TRGS methods stop. Figure 5 shows the plot of RSE versus IT (left) and RSE versus CPU (right) of Algorithm 2 applied to solve (1.1) with the sparse coefficient matrix $A$ named "divorce".

It can be seen from Table 8 that for sparse matrices listed in Table 7, TRGS achieves fast convergence with less IT and CPU(s) than that of RK, RGS, GTRK and RGS2 do. Similar results are shown in Table 9. Furthermore, from Figure 4, TRGS reaches the stop criteria with less iterations (left) and CPU time (right) than other methods for the "divorce" consistent sparse least squares problem (1.1). Similar results also appear in Figure 5.

Example 4.5 This example uses Algorithm 2 to restore a computer tomography (CT) image. We use the MATLAB function $paralleltomo(N, \theta, p)$ from Algebraic Iterative Reconstruction (ART) package in ^[31] to generate a large sparse matrix $A$ and the exact solution $x_*$, where $N = 35$, $\theta = 0^{\circ}:1.5^{\circ}:178^{\circ}$ and $p = 50$, then the size of $A$ is $5950\times 1225$ and the condition number $Cond(A) = 352.32$. We compute $\hat{b}$ by $\hat{b} = Ax_*$ and $b = \hat{b}+r$, where the noise $r$ is from the null space of the coefficient matrix $A^T$, i.e., $A^Tr = 0$. We set ordinary Gaussian white noise with noise levels $\delta = 0.01$ and the maximum iterative number is $5*10^6$. The TRGS is used to recover $x_*$ from $b$ and compared with the REK, RGS, PREK, RGS2 and TREK methods.

Figure 6 shows the recovered images by REK, RGS, PREK, RGS2, TREK and TRGS together with the original image and noised image with $\delta = 0.01$. Figure 7 shows the convergence of RSE versus IT (left) and RSE versus CPU(s) (right) of TRGS compared with other methods when $\delta = 0.01$.

It is shown from Figure 6 that all methods obtain well restored image, and Figure 7 shows that TRGS converges much faster than REK, RGS, PREK, RGS2 and TREK do when $\delta = 0.01$.

Example 4.6 In this example, we use Algorithm 2 to solve the famous Phillips ill-posed problem in ^[32], which comes from the Fredholm integral equation of first kind

on the square $[-6, 6]\times[-6, 6]$, where the kernel function is presented by $K(s, t) = \phi(s-t)$ with

and the right-hand side

The above problem is discretized to obtain the linear systems (1.1), where the coefficient matrix $A\in\mathbb{R}^{n\times n}$, exact solution vector $x_*\in\mathbb{R}^n$ and column vector $b\in\mathbb{R}^n$ are all generated by MATLAB function $phillips(n)$ in ^[32]. As we know, if the system matrix $A$ of (1.1) is ill-conditioned, one way to solve this problem is to use the Tikhonov regularization, that is

We set $n = 1024$, $\lambda = 0.01$ and Gaussian white noise level denoted by $\delta = \|r\|_2/||b||_2 = 1\%$ to obtain $b = Ax_*+r$. Here the condition number $Cond(A) = 4.1618e+10$ and the rank is 1024, which means (1.1) is a several linear ill-posed problem. It is worth noting that the system matrix satisfies the conditions of Lemma 1, Theorems 1 and 2, so the above methods will converge. In Figure 8, the (a) displays the approximation solution derived by RGS, RGS2 and TRGS together with the exact solution for the Phillips test problem when $\delta = 0.01$. The (b) and (c) show the convergence of RSE versus IT and RSE versus CPU(s) of TRGS compared with the RGS and RGS2 methods when $\delta = 0.01$.

We can see from Figure 8 that all the recovered solution by RGS, RGS2 and TRGS are close to the exact solution in (a), (b) and (c) show that TRGS converges much faster than the other two methods.

5.   Conclusions

A two-step randomized Guass-Seidel (TRGS) method for solving the linear least squares problems with tall and narrow coefficient matrix was presented. And the convergence analysis is provided when the coefficient matrix of (1.1) is of full column rank. This method does not need any columns pavings. Numerical examples for different cases show the superiority of the current method (TRGS) in this paper.

Acknowledgments

This work was supported by the Project of Department of Science and Technology of Sichuan Provincial (No. 2021ZYD0005), the Opening Project of Key Laboratory of Higher Education of Sichuan Province for Enterprise Informationalization and Internet of Things (No. 2020WZJ01), the Opening Project of Bridge Non-destruction Detecting and Engineering Computing Key Laboratory of Sichuan (Nos. 2021QYJ07, 2020QZJ03), the Scientific Research Project of Sichuan University of Science and Engineering (No. 2020RC24), and the Graduate student Innovation Fund (No. y2021101).

Conflict of interest

The authors declare no conflicts of interest regarding the publication of this paper.

Appendix

Proof of Theorem 2. By the definition of TRGS, one can obtain that

Then,

The following will estimate $\frac{A_{j_{k_1}}^T(b-Ay_k)}{\|A_{j_{k_1}}\|_2}$ and $\frac{A_{j_{k_2}}^T(b-Ay_k)}{\|A_{j_{k_2}}\|_2}$,

(i)

then,

(ii)

Substitute (A2) and (A3) into (A1), then,

Subtract $Ax_*$ from both sides of the equation, one can get

Since

and

then,

Let $\hat{y}_k$ and $\hat{x}_{k+1}$ be the first and second iterative solutions of $x_k$ obtained by single-step continuous execution of RGS2 algorithm, respectively. i.e.,

According to Theorem 1, one has

the last equation is obtained from (A4), then

Due to the orthogonality, one can get that,

So,

Next, we estimate the two conditional expectations on the right side of (A5).

(I) The first expectation

According to Theorem 1,

(II) The second expectation

By Lemma 3, $\|u_k\|^2_2-1 = -\mu_k^2$. Then,

By Lemma 4, it is established as follows

For any $\alpha, \beta, \theta, \eta\in\mathbb{R}$, we note the fact that

Then,

By Lemma 4, $(\lvert\alpha_{s, t}\lvert-\lvert\beta_{s, t}\lvert)^2\geq\gamma$. Then

i.e.,

Substitute (A6) and (A7) into (A5), one can obtain that

Thus,