Greedy Randomized Gauss-Seidel Method with Oblique Direction

Weifeng Li, Pingping Zhang

School of Science, Chongqing University of Posts and Telecommunications, Chongqing, China.

DOI: 10.4236/jamp.2023.114068                

Abstract

For the linear least squares problem with coefficient matrix columns being highly correlated, we develop a greedy randomized Gauss-Seidel method with oblique direction. Then the corresponding convergence result is deduced. Numerical examples demonstrate that our proposed method is superior to the greedy randomized Gauss-Seidel method and the randomized Gauss-Seidel method with oblique direction.

Keywords

Oblique Direction, Linear Least Squares Problem, Gauss-Seidel Method

1. Introduction

We focus on the linear least squares problem

$\underset{\beta \in {ℝ}^n}{\min }{\Vert y-X\beta \Vert }_2^2\mathrm{,}$ (1)

where $X\in {ℝ}^{m\times n}$ $\left(m\ge n\right)$ is of full column rank, $y\in {ℝ}^m$ and ${\Vert \cdot \Vert }_2$ denotes the Euclidean norm. There is a wide range of applications for the least squares problem in many fields, such as signal processing, image restoration, and so on. As we know, the coordinate descent method is an effective iteration method for (1). It applies the Gauss-Seidel method to the following equivalent normal Equation

$X^{\text{T}}X\beta =X^{\text{T}}y\mathrm{,}$ (2)

leading to the following formula

${\beta }_{k+1}={\beta }_k+\frac{X_{j_k}^{\text{T}}\left(y-X{\beta }_k\right)}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\mathrm{,}{j}_k=\left(k\mathrm{mod}n\right)+\mathrm{1,}$

where $X_{j_k}$ denotes the j_kth column of $X$ , $e_{j_k}$ is the j_kth unit coordinate vector and the superscript T denotes the transpose. The convergence of the Gauss-Seidel method highly depends on the order of the column selected in each step.

Inspired by the fancy work of Strohmer and Vershynin [1] , Leventhal and Lewis [2] proposed the randomized Gauss-Seidel (RGS) method and proved that it has an expected linear convergence rate. As Bai and Wu [3] pointed out an obvious flaw of the RGS method that the probability for selecting column will be a uniform column sampling if the coefficient matrix is scaled with a diagonal matrix. To tackle this problem, they proposed the greedy randomized coordinate descent (GRCD) or called the greedy randomized Gauss-Seidel (GRGS) method by adopting an effective probability for selecting the working column to capture larger entries of the residual vector with respect to (2). They showed that the GRGS method is significantly superior to the RGS method in terms of both theoretical analysis and numerical experiments. In addition, there are tons of attention about the Gauss-Seidel type methods [4] [5] [6] [7] . However, the convergence rate of the RGS method will be significantly reduced when the coefficient matrix columns are highly correlated. To improve the convergence rate, Wang et al. [8] proposed the randomized Gauss-Seidel method with oblique direction (RGSO) by combining two successive selected unit coordinate directions as the search direction

$d_k=e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\mathrm{.}$ (3)

They showed that, in terms of both theory and experiments, the RGSO method outperforms the RGS method. For more discussions about the oblique projection, we refer the readers to other literatures [9] [10] and their references.

However, the RGSO method still exists in the same flaw of the RGS method that the probability for selecting column will be a uniform column sampling if the coefficient matrix is scaled with a diagonal matrix. In addition, it is worth noting that the convergence rate of the GRGS method would significantly decrease when the coefficient matrix columns are close to linear correlation. For the above mentioned limitations, we present a greedy randomized Gauss-Seidel method with oblique direction (GRGSO) for solving (1), by combining the oblique direction and the GRGS method. In theory, it is proved that the iterative solutions generated by the GRGSO method can converge to the least squares solution ${\beta }_{\mathrm{<mo>\; *\; </mo>}}$ when the coefficient matrix is of full column rank. Numerical results show that compared with the RGSO and GRGS methods, the GRGSO method has a significant advantage over the iteration steps and computing time, especially when the coefficient matrix columns are highly correlated.

The organization of this paper is as follows. In Section 2, some notation and lemmas are introduced. In Section 3, we propose the GRGSO method for solving (1) and give its convergence analysis. Some examples are used to demonstrate the competitiveness of our proposed method in Section 4. Finally, we draw some brief conclusions in Section 5.

2. Notion and Preliminaries

In the beginning of this section, we give some notation. For a Hermitian positive definite matrix $B$ and a column vector $\beta$ with appropriate dimension, we denote ${\Vert \beta \Vert }_B^2=\langle B\beta \mathrm{,}\beta \rangle ={\Vert B^{\frac{1}{2}}\beta \Vert }_2^2$ and ${\beta }^{\left(i\right)}$ the ith entry of $\beta$ . For a given matrix $S\in {ℝ}^{m\times n}$ , ${\sigma }_{\min }\left(S\right)$ and ${\Vert S\Vert }_F$ denote the smallest nonzero singular value and the Frobenius norm of matrix $S$ , respectively. Let $r_k=y-X{\beta }_k$ and $s_k=X^{\text{T}}{r}_k$ , then $s_k^{\left(j\right)}=X_j^{\text{T}}{r}_k$ represents jth entry of $s_k$ . ${\beta }_{\mathrm{<mo>\; *\; </mo>}}$ is the optimal solution of the corresponding problem. We indicate by ${\mathbb{E}}_k$ the expected value conditional on the first k iterations, that is,

${\mathbb{E}}_k[\cdot ]=\mathbb{E}\left[\cdot \mathrm{<mo>\; |\; </mo>}{j}_0\mathrm{,}{j}_1\mathrm{,}\cdots \mathrm{,}{j}_{k-1}\right]\mathrm{,}$

where $j_t,t=0,1,\cdots ,k-1$ is the column selected at the t-th iteration.

In the following, we give a basic lemma.

Lemma 1 (See Bai and Wu [11] ) If the vector $u$ is in the column space of $A^{\text{T}}$ , it holds

${\Vert Au\Vert }_2^2\ge {\sigma }_{\min }^2\left(A\right){\Vert u\Vert }_2^2\mathrm{.}$

3. GRGSO Method

In this section, we design the GRGSO method for solving (1), by combining the oblique direction with the GRGS method. The pseudo-code of GRGSO method is listed in Table 1. The difference between the RGSO method and GRGSO method is the selection strategy. The RGSO method utilizes the random selection strategy. Specially, the RGSO method selects $j_{k+1}$ th column with probability $\frac{{\Vert X_{j_{k+1}}\Vert }_2^2}{{\Vert X\Vert }_F^2}$ in the numerical experiments, which can be equivalent to the uniform sampling if the Euclidean norms of all the columns of the matrix $X$ are same; while the GRGSO method aims to grasp the larger entries of the residual vector at each iteration. Compared with the GRGS method, our proposed method considers the oblique projection, which is expected to have better convergence performance in some cases of coefficient matrix columns being highly correlated.

Remark Let $t_k=\arg \underset{1\le j\le n}{\max }\left\{\frac{{\left|s_k^{\left(j\right)}\right|}_2^2}{{\Vert X_j\Vert }_2^2}\right\}$ , which implies $t_k\in V_k$ . Therefore, for all iterative step k, the index set $V_k$ generated by the GRGSO method is nonempty.

Remark In the GRGSO method, it holds that

$\begin{array}{c}{s}_{k+1}=X^{\text{T}}{r}_{k+1}=X^{\text{T}}\left(y-X{\beta }_{k+1}\right)\\ =X^{\text{T}}\left(y-X\left({\beta }_k+{\eta }_k^{\left(j_k\right)}{w}_{j_k}\right)\right)\\ =X^{\text{T}}\left(y-X{\beta }_k\right)-{\eta }_k^{\left(j_k\right)}{X}^{\text{T}}X{w}_{j_k}\\ =s_k-{\eta }_k^{\left(j_k\right)}{X}^{\text{T}}X{w}_{j_k}\mathrm{,}\left(k\ge 1\right)\end{array}$ (4)

Table 1. GRGSO method.

and

$\begin{array}{c}{s}_1=X^{\text{T}}{r}_1=X^{\text{T}}\left(y-X{\beta }_1\right)\\ =X^{\text{T}}\left(y-X\left({\beta }_0+\frac{s_0^{\left(j_1\right)}}{{\Vert X_{j_1}\Vert }_2^2}{e}_{j_1}\right)\right)\\ =X^{\text{T}}\left(y-X{\beta }_0\right)-\frac{s_0^{\left(j_1\right)}}{{\Vert X_{j_1}\Vert }_2^2}{X}^{\text{T}}{X}_{j_1}\\ =s_0-\frac{s_0^{\left(j_1\right)}}{{\Vert X_{j_1}\Vert }_2^2}{X}^{\text{T}}{X}_{j_1}\mathrm{.}\end{array}$ (5)

Therefore, the GRGSO method can be executed more effectively if $X^{\text{T}}X$ can be computed in an economical manner at the beginning.

Next, we give some lemmas which are useful to analyze the convergence of the GRGSO method.

Lemma 2 For the GRGSO method, we have

$s_k^{\left(j_k\right)}=0,\left(\forall k>0\right),$ (6)

$s_k^{\left(j_{k-1}\right)}=0,\left(\forall k>1\right).$ (7)

Proof. For $k=1$ , one has

$s_1^{\left(j_1\right)}=X_{j_1}^{\text{T}}\left(y-X{\beta }_1\right)=X_{j_1}^{\text{T}}y-X_{j_1}^{\text{T}}X\left({\beta }_0+\frac{s_0^{\left(j_1\right)}}{{\Vert X_{j_1}\Vert }_2^2}{e}_{j_1}\right)=0.$ (8)

For $k>1$ , we have

$\begin{array}{c}{s}_k^{\left(j_k\right)}=X_{j_k}^{\text{T}}\left(y-X{\beta }_k\right)=X_{j_k}^{\text{T}}y-X_{j_k}^{\text{T}}X\left({\beta }_{k-1}+{\eta }_{k-1}^{\left(j_{k-1}\right)}{w}_{j_{k-1}}\right)\\ =X_{j_k}^{\text{T}}y-X_{j_k}^{\text{T}}X{\beta }_{k-1}-{\eta }_{k-1}^{\left(j_{k-1}\right)}\left(X_{j_k}^{T}X{w}_{j_{k-1}}\right)={\tilde{s}}_{k-1}^{\left(j_k\right)}-\frac{{\tilde{s}}_{k-1}^{\left(j_k\right)}}{h_{j_{k-1}}}{h}_{j_{k-1}}\\ =\mathrm{0,}\end{array}$

which together with (8) proves (6).

By (6), it follows for $k>0$ that

$\begin{array}{c}{s}_{k+1}^{\left(j_k\right)}\mathrm{<mo>\; =\; </mo>}{X}_{j_k}^{\text{T}}\left(y-X{\beta }_{k+1}\right)=X_{j_k}^{\text{T}}y-X_{j_k}^{\text{T}}X\left({\beta }_k+{\eta }_k^{\left(j_k\right)}{w}_{j_k}\right)\\ =s_k^{\left(j_k\right)}-{\eta }_k^{\left(j_k\right)}\left(X_{j_k}^{\text{T}}X{w}_{j_k}\right)\\ =-{\eta }_k^{\left(j_k\right)}\left(X_{j_k}^{\text{T}}X\left(e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\right)\right)\\ =\mathrm{0,}\end{array}$

which leads to (7).

Remark From (6) and (7), it is obvious that in kth iteration, the GRGSO method dose not select $j_k$ , $j_{k-1}$ , which means $j_{k+1}\ne j_k\mathrm{,}{j}_{k-1}$ . Thus, the direction $w_{j_k}$ can be the combination of two unit coordinate directions. This is also the advantage of the GRGSO method compared with the RGSO method. Since the RGSO method randomly selects $j_{k+1}$ in the kth iteration, which can not avoid selecting $j_k$ and $j_{k-1}$ while the GRGSO method can avoid.

Lemma 3 For $h_{j_k}$ in the GRGSO method, it satisfies

$h_{j_k}={\Vert X{w}_{j_k}\Vert }_2^2={\Vert X_{j_{k+1}}\Vert }_2^2-\frac{{\left|X_{j_k}^{\text{T}}{X}_{j_{k+1}}\right|}^2}{{\Vert X_{j_k}\Vert }_2^2}\le \Delta {\Vert X_{j_{k+1}}\Vert }_2^2\mathrm{,}$

where $\Delta =\underset{i\ne j}{\max }{\sin }^2\langle X_i\mathrm{,}{X}_j\rangle$ .

Proof. Since $j_k\ne j_{k+1}$ , it holds that

$\begin{array}{c}{h}_{j_k}=X_{j_{k+1}}^{\text{T}}X\left(e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\right)={\Vert X_{j_{k+1}}\Vert }_2^2-\frac{{\left|X_{j_k}^{\text{T}}{X}_{j_{k+1}}\right|}^2}{{\Vert X_{j_k}\Vert }_2^2}\\ ={\Vert X_{j_{k+1}}\Vert }_2^2-\frac{{\Vert X_{j_k}\Vert }_2^2{\Vert X_{j_{k+1}}\Vert }_2^2{\cos }^2\langle X_{j_k}\mathrm{,}{X}_{j_{k+1}}\rangle }{{\Vert X_{j_k}\Vert }_2^2}\\ ={\sin }^2\langle X_{j_k}\mathrm{,}{X}_{j_{k+1}}\rangle {\Vert X_{j_{k+1}}\Vert }_2^2\le \Delta {\Vert X_{j_{k+1}}\Vert }_2^2\end{array}$

and

$\begin{array}{c}{\Vert X{w}_{j_k}\Vert }_2^2={\left(e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\right)}^{\text{T}}{X}^{\text{T}}X\left(e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\right)\\ ={\Vert X_{j_{k+1}}\Vert }_2^2-\frac{{\left|X_{j_k}^{\text{T}}{X}_{j_{k+1}}\right|}^2}{{\Vert X_{j_k}\Vert }_2^2}\mathrm{.}\end{array}$

Hence, we complete this proof.

Lemma 4 The iteration sequence ${\left\{{\beta }_k\right\}}_{k=0}^{\infty }$ generated by the GRGSO method satisfies

${\Vert X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\Vert X\left({\beta }_{k+1}-{\beta }_k\right)\Vert }_2^2\mathrm{,}k=\mathrm{0,1,2,}\cdots \mathrm{.}$

Proof. For $k=0$ , we have

$\begin{array}{c}{e}_{j_1}^{\text{T}}{X}^{\text{T}}X\left({\beta }_1-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)=X_{j_1}^{\text{T}}X\left({\beta }_0-{\beta }_{\mathrm{<mo>\; *\; </mo>}}+\frac{s_0^{\left(j_1\right)}}{{\Vert X_{j_1}\Vert }_2^2}{e}_{j_1}\right)\\ =X_{j_1}^{\text{T}}X{\beta }_0-X_{j_1}^{\text{T}}X{\beta }_{\mathrm{<mo>\; *\; </mo>}}+s_0^{\left(j_1\right)}\\ =0.\end{array}$

This means that the vector $X^{\text{T}}X\left({\beta }_1-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)$ is perpendicular to the vector $e_{j_1}$ . Since ${\beta }_1-{\beta }_0$ is parallel to $e_{j_1}$ , the vector $X^{\text{T}}X\left({\beta }_1-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)$ is perpendicular to ${\beta }_1-{\beta }_0$ .

For $k>0$ , It follows from Lemma 3 and Lemma 2 that

$\begin{array}{c}{w}_{j_k}^{\text{T}}{X}^{\text{T}}X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)=w_{j_k}^{\text{T}}{X}^{\text{T}}X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}+{\eta }_k^{\left(j_k\right)}{w}_{j_k}\right)\\ =w_{j_k}^{\text{T}}\left(-s_k+{\eta }_k^{\left(j_k\right)}{X}^{\text{T}}X{w}_{j_k}\right)\\ =-w_{j_k}^{\text{T}}{s}_k+{\eta }_k^{\left(j_k\right)}{\Vert X{w}_{j_k}\Vert }_2^2\\ =-{\left(e_{j_{k+1}}-\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{e}_{j_k}\right)}^{\text{T}}{s}_k+{\tilde{s}}_k^{\left(j_{k+1}\right)}\\ =-s_k^{\left(j_{k+1}\right)}+\frac{X_{j_k}^{\text{T}}{X}_{j_{k+1}}}{{\Vert X_{j_k}\Vert }_2^2}{s}_k^{\left(j_k\right)}+{\tilde{s}}_k^{\left(j_{k+1}\right)}=\mathrm{0,}\end{array}$

which means that the vector $X^{\text{T}}X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)$ is perpendicular to the vector $w_{j_k}$ . Since ${\beta }_{k+1}-{\beta }_k$ is parallel to $w_{j_k}$ , the vector $X^{\text{T}}X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)$ is perpendicular to ${\beta }_{k+1}-{\beta }_k$ . For all $k=\mathrm{0,1,}\cdots$ , it follows that

$\langle X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{,}X\left({\beta }_{k+1}-{\beta }_k\right)\rangle =\langle X^{\text{T}}X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\mathrm{,}{\beta }_{k+1}-{\beta }_k\rangle =\mathrm{0,}$

which together with Pythagoras theorem leads to the desired result.

Next, the convergence theory of the GRGSO method is deduced.

Theorem 5 For the least squares problem (1), the iteration sequence ${\left\{{\beta }_k\right\}}_{k=0}^{\infty }$ generated by the GRGSO method from any initial guess ${\beta }_0$ satisfies

$\mathbb{E}{\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\le {\displaystyle \underset{t=0}{\overset{k-1}{\prod }}}{\zeta }_t{\Vert X\left({\beta }_0-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\mathrm{,}\text{for}k=\mathrm{1,2,}\cdots \mathrm{,}$ (9)

where ${\zeta }_0=1-\frac{{\sigma }_{\min }^2\left(X\right)}{{\Vert X\Vert }_F^2}$ , ${\zeta }_1=1-\frac{1}{2}\left(\frac{1}{{\gamma }_1}+\frac{1}{{\Vert X\Vert }_F^2}\right)\frac{{\sigma }_{\min }^2\left(X\right)}{\Delta }$ , ${\zeta }_t=1-\frac{1}{2}\left(\frac{1}{{\gamma }_2}+\frac{1}{{\Vert X\Vert }_F^2}\right)\frac{{\sigma }_{\min }^2\left(X\right)}{\Delta }$ ( $t>1$ ). Here ${\gamma }_1=\underset{1\le i\le n}{\max }{\displaystyle \underset{\begin{array}{l}j=1\\ j\ne i\end{array}}{\overset{n}{\sum }}}{\Vert X_j\Vert }_2^2$ , ${\gamma }_2=\underset{\begin{array}{c}1\le i\mathrm{,}j\le n\\ i\ne j\end{array}}{\max }{\displaystyle \underset{\begin{array}{c}t=1\\ t\ne i\mathrm{,}j\end{array}}{\overset{n}{\sum }}}{\Vert X_t\Vert }_2^2$ and $\Delta$ is given as in Lemma 3.

Proof. By Lemma 2, it follows for $k>1$ that

$\begin{array}{c}{\delta }_k{\Vert X\Vert }_F^2=\frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\Vert s_k\Vert }_2^2}\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_k^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}+1\right)=\frac{1}{2}\left(\frac{\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_k^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}}\frac{{\Vert X_j\Vert }_2^2}{{\Vert X\Vert }_F^2}\frac{{\left|s_k^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}}+1\right)\\ =\frac{1}{2}\left(\frac{\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_k^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}}{{\displaystyle \underset{\begin{array}{c}j=1\\ j\ne j_k\mathrm{,}{j}_{k-1}\end{array}}{\overset{n}{\sum }}}\frac{{\Vert X_j\Vert }_2^2}{{\Vert X\Vert }_F^2}\frac{{\left|s_k^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}}+1\right)\ge \frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\displaystyle \underset{\begin{array}{c}j=1\\ j\ne j_k\mathrm{,}{j}_{k-1}\end{array}}{\overset{n}{\sum }}}{\Vert X_j\Vert }_2^2}+1\right)\\ \ge \frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\gamma }_2}+1\right)\mathrm{.}\end{array}$ (10)

For $k=1$ , it follows from (6) that

$\begin{array}{c}{\delta }_1{\Vert X\Vert }_F^2=\frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\Vert s_1\Vert }_2^2}\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_1^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}+1\right)=\frac{1}{2}\left(\frac{\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_1^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}}{{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{n}{\sum }}}\frac{{\Vert X_j\Vert }_2^2}{{\Vert X\Vert }_F^2}\frac{{\left|s_1^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}}+1\right)\\ =\frac{1}{2}\left(\frac{\underset{1\le j\le n}{\max }\left\{\frac{{\left|s_1^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}\right\}}{{\displaystyle \underset{\begin{array}{l}j=1\\ j\ne j_1\end{array}}{\overset{n}{\sum }}}\frac{{\Vert X_j\Vert }_2^2}{{\Vert X\Vert }_F^2}\frac{{\left|s_1^{\left(j\right)}\right|}^2}{{\Vert X_j\Vert }_2^2}}+1\right)\ge \frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\displaystyle \underset{\begin{array}{l}j=1\\ j\ne j_1\end{array}}{\overset{n}{\sum }}}{\Vert X_j\Vert }_2^2}+1\right)\\ \ge \frac{1}{2}\left(\frac{{\Vert X\Vert }_F^2}{{\gamma }_1}+1\right)\mathrm{.}\end{array}$ (11)

By Lemma 4, Lemma 3 and Lemma 1, for $k\ge 1$ , we have

$\begin{array}{l}{\mathbb{E}}_k{\Vert X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\\ ={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\mathbb{E}}_k{\Vert X\left({\beta }_{k+1}-{\beta }_k\right)\Vert }_2^2\\ ={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{\Vert {\eta }_k^{\left(j_k\right)}\left(X{w}_{j_k}\right)\Vert }^2\\ ={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2{\Vert X{w}_{j_k}\Vert }_2^2}{{\left|h_{j_k}\right|}^2}\\ ={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{h_{j_k}}\end{array}$

$\begin{array}{l}\le {\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{{\displaystyle \underset{j_{k+1}\in V_k}{\sum }}{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}\frac{{\left|{\tilde{s}}_k^{\left(j_{k+1}\right)}\right|}^2}{\Delta {\Vert X_{j_{k+1}}\Vert }_2^2}\\ \le {\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-\frac{{\delta }_k}{\Delta }{\Vert s_k\Vert }_2^2\\ ={\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2-\frac{{\delta }_k}{\Delta }{\Vert X^{\text{T}}X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\\ \le \left(1-{\delta }_k\frac{{\sigma }_{\min }^2\left(X\right)}{\Delta }\right){\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\mathrm{,}\end{array}$

which together with (11) and (10) can lead to

${\mathbb{E}}_1{\Vert X\left({\beta }_2-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\le {\zeta }_1{\Vert X\left({\beta }_1-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2$ (12)

and

${\mathbb{E}}_k{\Vert X\left({\beta }_{k+1}-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\le \zeta {\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\mathrm{,}k>\mathrm{1,}$ (13)

respectively. For $k=0$ , we can similarly get

$\mathbb{E}{\Vert X\left({\beta }_1-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\le {\zeta }_0{\Vert X\left({\beta }_0-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2\mathrm{.}$

Then taking the full expectation on both sides of (12) and (13) and by induction on the iteration index k, we can easily obtain (9). Thus, we complete the proof.

Remark Suppose that the upper bound for convergence rate of the GRGS method [3] is defined as

${\rho }_{GRCD}=1-\frac{1}{2}\left(\frac{1}{{\gamma }_1}+\frac{1}{{\Vert X\Vert }_F^2}\right){\sigma }_{\min }^2\left(X\right)\mathrm{.}$

Since $\Delta =\underset{i\ne j}{\max }{\sin }^2\langle X_i\mathrm{,}{X}_j\rangle \in \left(\mathrm{0,1}\right]$ , ${\gamma }_2<{\gamma }_1<{\Vert X\Vert }_F^2$ and $0<{\sigma }_{\min }^2\left(X\right)\le {\Vert X\Vert }_F^2$ , we have

$\zeta <{\zeta }_1\le {\rho }_{GRCD}<{\zeta }_0<1.$

This implies that the upper bound on the convergence factor of the GRGSO method is smaller, uniformly with respect to the iterative step k, than that of the GRCD method.

4. Numerical Experiments

Some examples are designed in this section to verify the effectiveness of the GRGSO method. Specially, the GRGSO method is compared with the GRGS method [3] , RGS method [2] and RGSO method [8] . We list the tested results of these methods in terms of the number of iteration steps (denoted by “IT”) and the running time in seconds (denoted by “CPU”).

The coefficient matrix $X\in {ℝ}^{m\times n}\left(m\ge n\right)$ in numerical experiments is from two sources. One is the random matrix whose entries are randomly taken from the interval $\left[c\mathrm{,1}\right]$ $\left(c\ge 0\right)$ by using the function rand in MATLAB. c being close to 1 implies that the matrix columns are highly correlated. Another is sparse matrices from the literature [12] listed in Table 2. We use $\text{cond}\left(X\right)$ to represent the condition number for $X$ , and define the density as

$\text{density}=\frac{\text{number}\text{of}\text{nonzero}\text{entries}\text{of}\text{an}m\times n\text{matrix}}{m\times n}\mathrm{.}$

We randomly generate the true vector ${\beta }_{\mathrm{<mo>\; *\; </mo>}}$ by utilizing the MATLAB function randn and construct the vector y by $y=X{\beta }_{\mathrm{<mo>\; *\; </mo>}}$ when the linear system is consistent; while for the inconsistent linear systems the right-hand side is set by $y=X{\beta }_{\mathrm{<mo>\; *\; </mo>}}+noise$ , where $noise\in null\left(X^{\text{T}}\right)$ . We take zero vector for the initial approximation of each iteration process. Since

${\Vert X\left({\beta }_k-{\beta }_{\mathrm{<mo>\; *\; </mo>}}\right)\Vert }_2^2={\Vert X{\beta }_{\mathrm{<mo>\; *\; </mo>}}+noise-r_k-X{\beta }_{\mathrm{<mo>\; *\; </mo>}}\Vert }_2^2={\Vert noise-r_k\Vert }_2^2,$

which was also used in the work of Wang et al. [8] , we terminate the iteration process once

$RSE=\frac{{\Vert noise-r_k\Vert }_2}{{\Vert y\Vert }_2}<{10}^{-6}\mathrm{.}$

We set “-” in the numerical tables if the number of iteration steps exceeds 300,000. All the results are averages from 20 repetitions. All experiments were implemented by using MATLAB (R2021b) on a computer with 2.30 GHz central processing unit (Intel(R) Core(TM) i7-10875H CPU), 16 GB memory.

Table 2. Sparse matrix properties of realistic problems [12] .

For the randomly generated matrices, with $c=0$ , the numerical results for consistent and inconsistent linear systems are listed in Table 3 and Table 4, respectively. It is easy to observe from Table 3 and Table 4 that the GRGSO method significantly outperforms the RGS, GRGS and RGSO methods in terms of both IT and CPU.

In the following, we compare these methods for solving (1) when the randomly generated matrix is with different c. We list the numerical results for the consistent and inconsistent systems in Table 5 and Table 6, respectively. From Table 5 and Table 6, it is easy to observe that the IT and CPU of the RGS and GRGS methods increase significantly with c increasing closer to 1. When c increases to 0.8, the IT of the RGS method exceeds the maximal number of iteration steps. And when c increases to 0.9, the IT of the GRGS method exceeds the maximal number of iteration steps. For the different c values, the methods with the oblique direction outperform the methods without the oblique direction. In addition, compared with the other methods, the GRGSO method performs best in terms of both IT and CPU.

For the full-rank sparse matrices from literature [12] , the numerical results for the consistent and inconsistent linear systems are listed in Table 7 and Table 8, respectively. It can be seen that for the matrix abtaha1 the performance of the GRGSO method is similar to that of the GRGS method in terms of both IT and CPU, whereas for the other matrices the GRGSO method significantly performs better in terms of both IT and CPU than the other methods.

Table 3. The consistent system with $c=0$ : different m impacts on IT and CPU.

Table 4. The inconsistent system with $c=0$ : different m impacts on IT and CPU.

Table 5. The consistent system with $X\in {ℝ}^{1000\times 100}$ : different c impacts on IT and CPU.

Table 6. The inconsistent system with $X\in {ℝ}^{1000\times 100}$ : different c impacts on IT and CPU.

Table 7. The consistent system: IT and CPU time of test methods for different sparse matrices.

Table 8. The inconsistent system: IT and CPU time of test methods for different sparse matrices.

5. Conclusion

In this manuscript, we construct the GRGSO method for the linear least squares problem. We have established the convergence analyses of the GRGSO method. Numerical experiments show that the GRGSO method is superior to the RGS, GRGS and RGSO methods in terms of both IT and CPU, especially when the coefficient matrix columns are highly correlated. It is natural to generalize the GRGSO method by introducing a relaxation parameter in its probability criterion. However, the choice of the optimal relaxation parameter is difficult in theory up to now. How to find the optimal relaxation parameter in theory is worthy of further study.

Conflicts of Interest

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

References

[1]	Strohmer, T. and Vershynin, R. (2009) A Randomized Kaczmarz Algorithm with Exponential Convergence. Journal of Fourier Analysis and Applications, 15, 262-278. https://doi.org/10.1007/s00041-008-9030-4
[2]	Leventhal, D. and Lewis, A.S. (2010) Randomized Methods for Linear Constraints: Convergence Rates and Conditioning. Mathematics of Operations Research, 35, 641-654. https://doi.org/10.1287/moor.1100.0456
[3]	Bai, Z.Z. and Wu, W.T. (2019) On Greedy Randomized Coordinate Descent Methods for Solving Large Linear Least-Squares Problems. Numerical Linear Algebra with Applications, 26, e2237. https://doi.org/10.1002/nla.2237
[4]	Du, K. (2019) Tight Upper Bounds for the Convergence of the Randomized Extended Kaczmarz and Gauss-Seidel Algorithms. Numerical Linear Algebra with Applications, 26, e2233. https://doi.org/10.1002/nla.2233
[5]	Liu, Y., Jiang, X.L. and Gu, C.Q. (2021) On Maximum Residual Block and Two-Step Gauss-Seidel Algorithms for Linear Least-Squares Problems. Calcolo, 58, Article No. 13. https://doi.org/10.1007/s10092-021-00404-x
[6]	Zhang, J.H. and Guo, J.H. (2020) On Relaxed Greedy Randomized Coordinate Descent Methods for Solving Large Linear Least-Squares Problems. Applied Numerical Mathematics, 157, 372-384. https://doi.org/10.1016/j.apnum.2020.06.014
[7]	Niu, Y.Q. and Zheng, B. (2021) A New Randomzied Gauss-Seidel Method for Solving Linear Least-Squares Problems. Applied Mathematics Letters, 116, Article ID: 107057. https://doi.org/10.1016/j.aml.2021.107057
[8]	Wang, F., Li, W.G., Bao, W.D. and Lv, Z.L. (2021) Gauss-Seidel Method with Oblique Direction. Results in Applied Mathematics, 12, Article ID: 100180. https://doi.org/10.1016/j.rinam.2021.100180
[9]	Li, W.G., Wang, Q., Bao, W.D. and Xing, L.L. (2022) Kaczmarz Method with Oblique Projection. Results in Applied Mathematics, 16, Article ID: 100342. https://doi.org/10.1016/j.rinam.2022.100342
[10]	Wang, F., Li, W.G., Bao, W.D. and Liu, L. (2022) Greedy Randomized and Maximal Weighted Residual Kaczmarz Methods with Oblique Projection. Electronic Research Archive, 30, 1158-1186. https://doi.org/10.3934/era.2022062
[11]	Bai, Z.Z. and Wu, W.T. (2018) On Greedy Randomized Kaczmarz Method for Solving Large Sparse Linear Systems. SIAM Journal on Scientific Computing, 40, A592-A606. https://doi.org/10.1137/17M1137747
[12]	Davis, T.A. and Hu, Y. (2011) The University of Florida Sparse Matrix Collection. ACM Transactions on Mathematical Software, 38, 1-25. https://doi.org/10.1145/2049662.2049663