Randomized symmetric Gauss-Seidel method for solving linear least squares problems

1.   Introduction

We consider solving the linear least squares problems given by

where $A\in \mathbb{R}^{m\times n}$ has full column rank and $m\geq n$. Problem (1.1) is prevalent in several fields, including ridge regression, machine learning, optimal control, and others. To address this problem in a resourceful and efficient manner, a substantial body of research has been conducted on iterative methods ^[1,2,3]. Notably, the randomized Kaczmarz method ^[4] and the randomized Gauss-Seidel method ^[5,6,7,8,9] have amassed considerable interest due to their capacity to handle large volumes of data.

Strohmer and Vershynin ^[4] put forth a randomized version of the Kaczmarz method. This method uniquely selects a row in proportion to the squared Euclidean norm, resulting in a fast convergence. Inspired by this, Leventhal and Lewis ^[10] developed a randomized Gauss-Seidel (RGS) method that samples a column of $A$ based on an appropriately chosen probability. Subsequently, Ma et al. ^[11] established a comprehensive convergence theory.

To expedite convergence, researchers have extensively investigated the two-step Gauss-Seidel method. Liu et al. ^[12] proposed the 2SGS method, a deterministic iteration scheme based on the maximum residual rule. Liao et al. ^[13] introduced RGS2, a method that merges two single-step iterations into one, sampling two distinct indices simultaneously. In the same literature, the TRGS method, an enhanced two-step approach, projects the approximate solution onto the solution space using two random columns. Mustafa and Saha ^[14] devised D2RGS, a two-dimensional coordinate descent method employing uniform sampling to randomly select two distinct columns of the coefficient matrix in each iteration.

In this paper, we concentrate on the result found by Niu and Zheng ^[8], which presents a novel randomized Gauss-Seidel (NRGS) method as follows:

where $A_{i_{k}}$ is the $i_{k}$-th column of $A$, $e_{i_{k}}$ is the ${i_{k}}$-th column of the $n\times n$ identity matrix, $r_{k} = b-Ax_{k}$, and the column index is chosen with probability $p_{i_{k}} = \frac{|A_{i_{k}}^\top r_{k}|^2}{\|{A^\top r_{k}}\|_2^2}$. The convergence result was given by

where $\beta = \|{A}\|_F^2-\min_{ 1\leq i \leq n}\|{A_i}\|_2^2$, $A^{\dagger}$ is the Moore-Penrose pseudoinverse of $A$, and $A^{\dagger}b$ is the unique solution of (1.1).

For solving (1.1) more efficiently, we consider projecting the update residual to two random symmetric columns. Borrowing the idea of the probability criterion given by NRGS in ^[8], we introduce Random Symmetric Gauss-Seidel (RSGS) method.

The organization of this paper is as follows. In Section 2, we present the related work. The RSGS method and its convergence analysis are presented in Section 3. Several numerical examples are displayed in Section 4 to show that our proposed RSGS method performs better than NRGS and other methods compared. Finally, the conclusion is drawn in Section 5.

2.   Related work

A large number of methods have been designed to solve linear least squares problems. Here, we review only the most relevant and recent research works.

Niu and Zheng ^[8] proposed a single-step Gauss-Seidel method, namely NRGS. This algorithm adopts an effective probabilistic criterion, allowing it to efficiently capture the larger elements within the residual vector. Inspired by its probabilistic criterion, we extend this algorithm to the two-dimensional scenario. While ensuring the simplification of the sampling method, we achieve faster convergence rates.

To accelerate the convergence, numerous researchers have conducted studies on the two-step Gauss-Seidel method. Liu et al. ^[12] proposed a two-step iteration Gauss-Seidel deterministic method named the 2SGS method, which is based on the maximum residual rule. As a deterministic algorithm, 2SGS is particularly suitable for medium-sized problems. However, for large-scale datasets, it is necessary to investigate stochastic algorithms. 2SGS is a highly effective deterministic approach, serving as an inspiring foundation for our subsequent design of efficient stochastic algorithms in large-scale datasets.

Liao et al. ^[13] introduced the two-step iteration method RGS2, which effectively combines two single-step iterations into one step, where two different indices are sampled simultaneously. The advantage of this approach is the avoidance of the possibility of repeating the same index in two consecutive single-step iterations, thereby improving convergence efficiency. In the same literature, the TRGS is the improved two-step method that projects the approximate solution onto the solution space by two random columns. Differently, each iteration of our method employs a stochastic two-dimensional coordinate space least squares approach. Intuitively, under the same conditions, the two-dimensional minimization method outperforms two-step minimization approach (i.e., the combination of two separated singe steps).

Mustafa and Saha ^[14] developed a two-dimensional coordinate descent method D2RGS, which employs uniform sampling to randomly select two different columns of the coefficient matrix in each iteration. While two-dimensional uniform sampling can save time compared to non-uniform sampling, it fails to consider the importance of different columns. The challenge lies in balancing the fast convergence of the algorithm with the efficiency of sampling. We aim to achieve local minimization in the two-dimensional coordinate space through symmetric sampling. Our method requires only non-uniform sampling of one index at a time while obtaining two indices. This approach is both simple to implement and efficient.

3.   The randomized symmetric Gauss-Seidel method

In this section, we introduce the randomized symmetric Gauss-Seidel method (RSGS), which can also be explained as the randomized symmetric coordinate projection method. The iteration scheme is given by

where $\alpha_k$ and $\beta_k$ are parameters which are chosen dynamically such that

where $r_{k+1} = b-Ax_{k+1}$. Substituting (3.1) into (3.2), we have

Thus,

and

Here, $A_j$ and $e_j$ are the $j$-th column of $A$, and the identity matrix $I_n$, respectively. The sampled index $i_k\in [n]\triangleq \{1, 2, \cdots, n\}$ is chosen with the probability

Based on this construction and criterion, we give Algorithm 1.

To discuss the convergence, we first give two Lemmas.

Lemma 1. Let $\alpha_k, \beta_k$ satisfy (3.4), (3.5) and (3.3). Then it holds that

and

where if $i_k = n-i_k+1$, equalities hold in (3.7) and (3.8).

Proof. Substituting the first equality of (3.3) into the right-hand side of (3.7), we obtain

According to Cauchy-Schwarz inequality $A_{i_k}^\top A_{n-i_k+1}\leq \|A_{i_k}\|_2\|A_{n-i_k+1}\|_2$, (3.7) is proved.

Similarly, substituting the second inequality of (3.3) into the right-hand side of (3.8) and by Cauchy-Schwarz inequality, (3.8) is proved.

Note that if $i_k = n-i_k+1$, then $A_{i_k} = A_{n-i_k+1}$ and $A_{i_k}^\top A_{n-i_k+1} = \|A_{i_k}\|_2\|A_{n-i_k+1}\|_2$, then equalities hold in (3.7) and (3.8).

The proof is finished. □

Lemma 2. ^[15] The following inequality holds

for $x_i\geq 0, y_i > 0, l > 0, i = 1, 2, \ldots, n.$ The equality holds if and only if $\frac{x_1}{y_1} = \cdots = \frac{x_n}{y_n}$.

The convergence of Algorithm 1 is given in the following theorem.

Theorem 1. Let $A\in \mathbb{R}^{m\times n}$ be of full column rank, and $m\geq n$. For the linear least squares problems (1.1), the iteration series $\{x_{k}\}_{k = 0}^{\infty}$ generated by Algorithm 1 converges in expectation to the unique solution $x_* = A^{\dagger}b$. Furthermore, we have the following results:

● If $n$ is even, for $k\geq 0$, we have

● If $n$ is odd, we have

where $\gamma_1 = \|{A}\|_{F}^2-\min_{1\leq i\leq n}\|{A_i}\|_2^2$, and $\gamma_2 = \|{A}\|_{F}^2-2\min_{1\leq i\leq n}\|{A_i}\|_2^2$.

Proof. According to Algorithm 1, the $k+1$-th $(k\geq 0)$ iteration can be represented as the optimization problem:

which is equivalent to the following orthogonal projection,

Then, according to (3.3)–(3.5), for $k\geq 0$, we have

where the second inequality is due to the fact that, for $k\geq 0$, $|A_{i_k}^\top A_{n-i_k+1}|\leq\|{A_{i_k}}\|_2\|{A_{n-i_k+1}}\|_2$. Taking the conditional expectation and according to (3.9), for $k\geq 0$, it holds that

Now, we estimate the lower bound of the second term in the above last equality. By (3.3)–(3.5), for $k\geq 1$, we have

where the first inequality is obtained by Lemma 1, the second inequality is due to Lemma 2 and $\text{span}\{A_{i_{k-1}}, A_{i_{n-i_{k-1}+1}}\} \perp r_k$. Thus,

where

Note that, if $n$ is odd and $k\geq 2$, we have

where the two inequalities next to the last are due to that $r_{k-1}\bot A_{\frac{n+1}{2}}$, while the last inequality is obtained by (3.11). If $n$ is even, then

By taking the full expectation on both sides of (3.10), for $k\geq 2$, we get that

Note that

where the last inequality is obtained by Theorem 2.2 in ^[8].

By induction on $k$, if $n$ is even, for $k\geq 0$, we have

if $n$ is odd, we have

the proof is completed. □

Remark 1. Since $\|{A}\|_F^2-2\min\{\|{A_i}\|_2^2, i\in[n]\} < \|{A}\|_F^2-\min\{\|{A_i}\|_2^2, i\in[n]\}$, our method is superior to NRGS in theory.

4.   Numerical examples

In this section, we present several examples that utilize two groups of real coefficient matrices. These examples are designed to compare the effectiveness of RSGS and NRGS, RGS2 ^[13], TRGS ^[13], and D2RGS ^[14] methods when solving Problem (1.1). The first group of matrices is chosen from the Florida sparse matrix collection ^[16] and is listed in Table 1, where for the matrices nemsafm, df2177 and bibd_16_8, we set $A$ as their transpose. The second group consists of randomly generated dense matrices, produced using MATLAB's randn function. The inconsistent linear system is constructed by setting $b = Ax+r$, where $x$ is a vector with entries generated from a standard normal distribution, and the residual $r$ belongs to the null space of $A^\top$, which is derived using the null function in MATLAB. For all computations, the initial point is set as $x_0 = 0$, and the stopping criterion is:

In order to evaluate and compare the performance of RSGS, NRGS and other baseline algorithms, we graph the relative error (err) against the metrics of IT and CPU times. These graphs are depicted in Figures 1 and 2, respectively. As can be observed from the figures, RSGS demonstrates more efficiency than the NRGS, RGS2 ^[13], TRGS ^[13], and D2RGS ^[14] methods.

5.   Conclusions

We proposed a randomized symmetric Gauss-Seidel method for solving linear least squares problems with the nonuniform sampling on the probability criterion (3.6). Our theoretical analysis indicates that RSGS converges when the coefficient matrix has full column rank. Furthermore, numerical experiments demonstrate that RSGS outperform the baseline algorithms.

Author contributions

Fan Sha: Write original draft, Study conception, Formal analysis; Jianbing Zhang: Supervision, Revision. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is supported by National Key R&D Program of China (Grant No. 2022YFA1004403) and by funding from oppo: Multi-modal information enhanced Image Caption.

Conflict of interest

All authors declare no conflicts of interest in this paper.