A two-dimensional randomized extended Gauss-Seidel algorithm for solving least squares problems

Abstract

We study a two-dimensional coordinate descent method to solve large linear least squares problems expanding on the method presented by Leventhal and Lewis. For an overdetermined system with full column rank, they proved its convergence in expectation, regardless of its consistency. In our work, we present a block version of the same. We also present an update on the extension done by Ma et al. to address non-full rank linear systems or underdetermined linear systems. Convergence is analyzed for the stated methods, and numerical experiments are provided to demonstrate their efficiency.

Appendix A   Convergence of R(E)BGS algorithm

Suppose, the REBGS method is implemented as per the notations used in Sect. 1.3. To obtain a $\left( \tilde{p}, \tilde{\alpha }, \tilde{\beta } \right)$ column paving of A, we randomly choose columns of A in such a way that $\tilde{p}$ number of selections will exhaust all the columns. A naive way may be to randomly choose $\dfrac{n}{\tilde{p}}$ columns for $\tilde{p}$ times without replacement, for some divisor $\tilde{p}$ of n. The row pavings can simply be seen as a column partition executed on $A^{T}$. The convergence results depend on the positive numbers $\gamma =\left( 1-\dfrac{\sigma _{min}^2A}{p\beta }\right)$, $\overline{\gamma } =\left( 1-\dfrac{\sigma _{min}^2A}{\tilde{p}\tilde{\beta }}\right)$. For a matrix $A\in \mathbb {R}^{m,n}$, with a $\left( p, \alpha , \beta \right)$ row paving $\rho$ and for any $z\in R(A^T)$, $\rho _i\in \rho$, we have [2, 14],

$$\begin{aligned} E\Vert \left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) z\Vert ^2\le \gamma \Vert z\Vert ^2. \end{aligned}$$

(A1)

Let $x_k$ be the $k^{th}$ iterate of the REBGS algorithm with arbitrary $x_0$. Then, in exact arithmetic,

$$\begin{aligned} E\Vert Ax_k-AA^{\dagger }b\Vert ^2\le \overline{\gamma }^k \Vert Ax_0-AA^{\dagger }b\Vert ^2. \end{aligned}$$

(A2)

We have

$$\begin{aligned} Ax_k-AA^{\dagger }b&=Ax_{k-1}-AA^{\dagger }b-A_{\tau _j}A_{\tau _j}^{\dagger }\left( Ax_{k-1}-b\right) \\&=\left( I-A_{\tau _j}A_{\tau _j}^{\dagger }\right) \left( Ax_{k-1}-AA^{\dagger }b\right) , \end{aligned}$$

where the second equality follows from the fact that $A_{\tau _j}A_{\tau _j}^{\dagger }b=A_{\tau _j}A_{\tau _j}^{\dagger }AA^{\dagger }b$.

Since, $A_{\tau _j}A_{\tau _j}^{\dagger }=\left( A_{\tau _j}^T\right) ^{\dagger }A_{\tau _j}^T$ and $Ax_{k-1}-AA^{\dagger }b\in R(A)$, we use (A1) to get

$$E_{k-1}\Vert Ax_k-AA^{\dagger }b\Vert ^2\le \overline{\gamma } \Vert Ax_{k-1}-AA^{\dagger }b\Vert ^2.$$

Using law of total expectation and then rolling out the recurrences, we obtain (A2), which in fact shows the convergence of the RBGS algorithm for full column rank systems.

The next result shows the convergence of the REBGS method which says that for $A\in \mathbb {R}^{m,n}$ and $b\in \mathbb {R}^{m}$, if $z_k$ denotes the $k^{th}$ iterate of the REBGS algorithm with arbitrary $x_0$ and $z_0 \in R(A^T)$, then, in exact arithmetic,

$$\begin{aligned} E\Vert z_k-A^{\dagger }b\Vert ^2\le \gamma ^k \Vert z_{0}-A^{\dagger }b\Vert ^2+r\left( \gamma ,\overline{\gamma }\right) \dfrac{\Vert Ax_0-AA^{\dagger }b\Vert ^2}{\alpha p} , \end{aligned}$$

(A3)

where for any $k\ge 1$,

$$r\left( \gamma ,\overline{\gamma }\right) = {\left\{ \begin{array}{ll} \overline{\gamma }\dfrac{\overline{\gamma }^k-\gamma ^k}{\overline{\gamma }-\gamma }, &{} \text { if } \overline{\gamma }\ne \gamma ,\\ k\gamma ^k, &{} \text { if }\overline{\gamma }= \gamma . \end{array}\right. } $$

We have

$$\begin{aligned} z_k-A^{\dagger }b&=z_{k-1}+A_{\rho _i}^{\dagger }\left( (Ax_k)_{\rho _i}-A_{\rho _i}z_{k-1}\right) -A^{\dagger }b\\&=\left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) z_{k-1} +A_{\rho _i}^{\dagger }A_{\rho _i}x_k-A^{\dagger }b\\&=\left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) \left( z_{k-1}-A^{\dagger }b\right) +A_{\rho _i}^{\dagger }A_{\rho _i}\left( x_k-A^{\dagger }b\right) . \end{aligned}$$

It follows from the orthogonality of the range spaces of $I-A_{\rho _i}^{\dagger }A_{\rho _i}$ and $A_{\rho _i}^{\dagger }A_{\rho _i}$, that

$$\begin{aligned} \Vert z_k-A^{\dagger }b\Vert ^2=\Vert \left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) \left( z_{k-1}-A^{\dagger }b\right) \Vert ^2+\Vert A_{\rho _i}^{\dagger }A_{\rho _i}\left( x_k-A^{\dagger }b\right) \Vert ^2. \end{aligned}$$

(A4)

As $z_0$, $A^{\dagger }b$ $\in R(A^T)$, it can be induced that $z_{k-1}-A^{\dagger }b\in R(A^T)$, and so by relation (A1), we have

$$E_k\Vert \left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) \left( z_{k-1}-A^{\dagger }b\right) \Vert ^2\le \gamma \Vert z_{k-1}-A^{\dagger }b\Vert ^2.$$

By law of total expectation,

$$\begin{aligned} E\Vert \left( I-A_{\rho _i}^{\dagger }A_{\rho _i}\right) \left( z_{k-1}-A^{\dagger }b\right) \Vert ^2\le \gamma E\Vert z_{k-1}-A^{\dagger }b\Vert ^2. \end{aligned}$$

(A5)

Again,

$$\begin{aligned} \Vert A_{\rho _i}^{\dagger }A_{\rho _i}\left( x_k-A^{\dagger }b\right) \Vert ^2&\le \sigma _{max}^2\left( A_{\rho _i}^{\dagger } \right) \Vert \left( Ax_k-AA^{\dagger }b\right) _{\rho _i}\Vert ^2\\&=\dfrac{ \Vert \left( Ax_k-AA^{\dagger }b\right) _{\rho _i}\Vert ^2}{\sigma _{min}^2\left( A_{\rho _i}^{\dagger } \right) }\\&\le \dfrac{\Vert \left( Ax_k-AA^{\dagger }b\right) _{\rho _i}\Vert ^2}{\alpha }, \end{aligned}$$

where the first inequality and the second equality result from properties of singular values and the third inequality follows from the paving relations. Therefore,

$$E_k\Vert A_{\rho _i}^{\dagger }A_{\rho _i}\left( x_k-A^{\dagger }b\right) \Vert ^2\le \dfrac{1}{\alpha p}\Vert Ax_k-AA^{\dagger }b\Vert ^2.$$

Taking expectation on both sides and using (A2), we get

$$\begin{aligned} E\Vert A_{\rho _i}^{\dagger }A_{\rho _i}\left( x_k-A^{\dagger }b\right) \Vert ^2\le \dfrac{\overline{\gamma }^k}{\alpha p} \Vert Ax_0-AA^{\dagger }b\Vert ^2. \end{aligned}$$

(A6)

It follows from (A4), (A5), and (A6) that

$$\begin{aligned} E\Vert z_k-A^{\dagger }b\Vert ^2&\le \gamma E\Vert z_{k-1}-A^{\dagger }b\Vert ^2+\dfrac{\overline{\gamma }^k}{\alpha p} \Vert Ax_0-AA^{\dagger }b\Vert ^2\\&\le \gamma ^2 E\Vert z_{k-2}-A^{\dagger }b\Vert ^2+\dfrac{\gamma \overline{\gamma }^{k-1}+\overline{\gamma }^{k}}{\alpha p} \Vert Ax_0-AA^{\dagger }b\Vert ^2\\&~~~~~~~~~\vdots ~~~~~~~~~\vdots ~~~~~~~~~\vdots \\&\le \gamma ^k E\Vert z_{0}-A^{\dagger }b\Vert ^2+\dfrac{1}{\alpha p}\sum _{i=0}^{k-1}\overline{\gamma }^{k-i}\gamma ^i\Vert Ax_0-AA^{\dagger }b\Vert ^2\\&=\gamma ^k \Vert z_{0}-A^{\dagger }b\Vert ^2+r\left( \gamma ,\overline{\gamma }\right) \dfrac{\Vert Ax_0-AA^{\dagger }b\Vert ^2}{\alpha p} . \end{aligned}$$

Hence, the (A3) holds.