A New Cyclic Gradient Method Adapted to Large-Scale Linear Systems

TL;DR

This paper introduces a novel cyclic gradient method tailored for large-scale linear systems, demonstrating finite termination in two dimensions and linear convergence in higher dimensions, with practical advantages shown through experiments.

Contribution

The paper presents a new gradient algorithm with proven finite termination in 2D and linear convergence in higher dimensions, outperforming existing methods in large-scale problems.

Findings

Finite termination in two dimensions.

R-linear convergence in any dimension.

Outperforms existing methods in large-scale experiments.

Abstract

This paper proposes a new gradient method to solve the large-scale problems. Theoretical analysis shows that the new method has finite termination property for two dimensions and converges R-linearly for any dimensions. Experimental results illustrate first the issue of parallel implementation. Then, the solution of a large-scale problem shows that the new method is better than the others, even competitive with the conjugate gradient method.

Tables

Table 1. Table 1: Gradient methods with different residual threshold

	${10}^{-1}$	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$
CG	$58$	$735$	$2617$	4251	5786	7535
CY	13	208	1153	4275	6181	\
CSD	$27$	$212$	1470	$4357$	$6713$	\
CBB	$31$	$241$	$1965$	$7534$	\	\
DY	15	200	$1595$	$6415$	\	\
BB1	$356$	$984$	$2494$	$4612$	$8633$	\
SD	$61$	$5773$	\	\	\	\

Equations

$$Ax=b,$$

$$f(x)=\frac{1}{2}{x}^{⊺}Ax-b^{⊺}x.$$

$$x_{k+1}=x_k-{\alpha }_k{g}_k,k=0,1,\dots ,$$

$${\alpha }_k^{\text{SD}}=\frac{g_k^{⊺}{g}_k}{g_k^{⊺}A{g}_k}.$$

$${\alpha }_k^{\text{BB1}}=\frac{g_{k-1}^{⊺}{g}_{k-1}}{g_{k-1}^{⊺}A{g}_{k-1}}.$$

$${\alpha }_k^{\text{BB2}}=\frac{g_{k-1}^{⊺}A{g}_{k-1}}{g_{k-1}^{⊺}{A}^2{g}_{k-1}}.$$

$${\alpha }_k^{\text{Y}}=\frac{2}{\sqrt{{\left(\frac{1}{{\alpha }_{k-1}^{\text{SD}}}-\frac{1}{{\alpha }_k^{\text{SD}}}\right)}^2+\frac{4g_k^{⊺}{g}_k}{s_{k-1}^{⊺}{s}_{k-1}}}+\frac{1}{{\alpha }_{k-1}^{\text{SD}}}+\frac{1}{{\alpha }_k^{\text{SD}}}},$$

$${\alpha }_k^{\text{DY}}=\{\begin{array}{ll}{\alpha }_k^{\text{SD}},& k\mathrm{mod}4<2,\\ {\alpha }_k^{\text{Y}},& \text{otherwise}.\end{array}$$

$$\bar{k}=\max \{0,k-m\},$$

$${\alpha }_k^{\text{GMR}}=\frac{g_{\tau (k)}^{⊺}{A}^{\rho (k)}{g}_{\tau (k)}}{g_{\tau (k)}^{⊺}{A}^{\rho (k)+1}{g}_{\tau (k)}},$$

$$\tau (k)\in \left\{\bar{k},\bar{k}+1,\dots ,k-1,k\right\},$$

$$\rho (k)\in \left\{q_1,\dots ,q_m\right\},q_j\ge 0,$$

$${\alpha }_k^{\text{CSD}}=\{\begin{array}{ll}{\alpha }_k^{\text{SD}},& k\mathrm{mod}m=0,\\ {\alpha }_{k-1},& \text{otherwise}.\end{array}$$

$${\alpha }_k^{\text{CBB}}=\{\begin{array}{ll}{\alpha }_k^{\text{BB1}},& k\mathrm{mod}m=0,\\ {\alpha }_{k-1},& \text{otherwise}.\end{array}$$

$${\alpha }_k^{\text{YB}}=\{\begin{array}{ll}{\alpha }_k^{\text{SD}},& k\mathrm{mod}3=0\text{ or }2,\\ {\alpha }_k^{\text{Y}},& k\mathrm{mod}3=1.\end{array}$$

$$\forall m\in \mathbb{N},\text{ if }k\mathrm{mod}(3+m)>2,\text{ then }{\alpha }_k={\alpha }_{k-1}.$$

$${\alpha }_k^{\text{CY}}=\{\begin{array}{ll}{\alpha }_k^{\text{Y}},& k\mathrm{mod}(l+m+2)=1,\\ {\alpha }_k^{\text{SD}},& k\mathrm{mod}(l+m+2)<l+2,\\ {\alpha }_{k-1},& \text{otherwise},\end{array}$$

$$A=\text{diag}({\lambda }_1,\dots ,{\lambda }_n),$$

$$1={\lambda }_1\le \cdots \le {\lambda }_n.$$

$$G(k,\mu )=\sum _{i=1}^{\mu }{g}_{i,k}^2,$$

$${\left(\frac{1}{{\alpha }_{k-1}^{\text{SD}}}+\frac{1}{{\alpha }_k^{\text{SD}}}\right)}^{-1}<{\alpha }_k^{\text{Y}}<\min \left\{{\alpha }_{k-1}^{\text{SD}},{\alpha }_k^{\text{SD}}\right\},$$

$${\lambda }_1\le \frac{1}{{\alpha }_k^{\text{SD}}}<\frac{1}{{\alpha }_k^{\text{Y}}}<\frac{1}{{\alpha }_{k-1}^{\text{SD}}}+\frac{1}{{\alpha }_k^{\text{SD}}}\le 2{\lambda }_n.$$

$$G(k,\mu )\le ϵ,g_{\mu +1,k}^2\ge M_2ϵ,$$

$$\begin{array}{rl}\frac{1}{{\alpha }_k^{\text{Y}}}& >\frac{1}{{\alpha }_k^{\text{SD}}}=\frac{g_k^{⊺}A{g}_k}{g_k^{⊺}{g}_k}=\frac{{\sum }_{i=1}^n{\lambda }_i{g}_{i,k}^2}{{\sum }_{i=1}^n{g}_{i,k}^2}\\ & \ge \frac{{\lambda }_{\mu +1}{\sum }_{i=\mu +1}^n{g}_{i,k}^2}{{\sum }_{i=1}^{\mu }{g}_{i,k}^2+{\sum }_{i=\mu +1}^n{g}_{i,k}^2}\\ & \ge \frac{{\lambda }_{\mu +1}}{\frac{{\sum }_{i=1}^{\mu }{g}_{i,k}^2}{g_{\mu +1,k}^2}+1}\\ & \ge \frac{{\lambda }_{\mu +1}}{\frac{ϵ}{2ϵ}+1}=\frac{2}{3}{\lambda }_{\mu +1}\end{array}$$

Full text

Keywords: gradient methods; linear systems; large-scale problems; Barzilai-Borwein methods

1 Introduction

We are interested in investigating new gradient methods for the solution of linear system

[TABLE]

where $A\in\mathbb{R}^{n\times n}$ is symmetric positive definite (SPD) and $b\in\mathbb{R}^{n}$. This problem is equivalent to the minimization of a convex quadratic function

[TABLE]

Gradient methods generate a sequence of the form

[TABLE]

where $g_{k}=Ax_{k}-b$. It is well known that the steepest descent (SD) method [3] performs poor in most cases, where the steplength can be written as follows

[TABLE]

The iterates generated tend to asymptotically alternate between two directions [1]. In contrast, the conjugate gradient (CG) method [11] is often the method of choice that will terminate in at most $n$ iterations. It is very attractive because of its high efficiency and low storage requirement. Nonetheless, CG iteration depends strongly on the search of direction calculation, i.e., any derivation such as round-off errors can seriously degrade performance [9].

In the past several decades, a renewed interest for gradient methods has appeared since Barzilai and Borwein [2] proposed two efficient nonmonotone steplengths

[TABLE]

[TABLE]

The motivation consists in approximating the Hessian and imposing some quasi-Newton properties. Some theories and experiments have shown that BB methods have good performance and are competitive with CG methods when low accuracy is required or small perturbation exists [9]. The convergence has been proven by Raydan [19]. Furthermore, Friedlander et al. [10] provided a general framework under the name of “gradient method with retards” that SD and BB both belong to it, as well as several alternate methods proposed later [4, 5, 21].

Motivated by the two-dimensional finite termination property, Yuan [20] provided a somewhat complicated steplength

[TABLE]

where $s_{k-1}=x_{k}-x_{k-1}$. He gave two algorithms and some variants were investigated further by Dai and Yuan [6]. Among these methods, the second variant (DY) is the most efficient one according to the experiments in [6], where iterates are generated of the form

[TABLE]

In this paper, we address the properties of cyclic gradient methods, especially their parallel behavior. We propose a new algorithm based on the Yuan steplength, which has also the two-dimensional finite termination property. In the next section, we introduce the cyclic gradient methods and propose our new steplength. In Section 3, we give the convergence results of the new method. Some numerical results are presented in Section 4. Finally, a concluding remark is shown in Section 5.

2 Cyclic Gradient Methods

Friedlander et al. [10] proposed an ingenious framework that gives rise to a great number of potentially efficient algorithms. Firstly, assume that $m\in\mathbb{N}$ represents retard that allows to employ the information from previous iterations. Let

[TABLE]

then a collection of possible choices of steplength can be set as follows

[TABLE]

where

[TABLE]

and

[TABLE]

where $k\in\mathbb{N}$. The next theorem summarizes the convergence result in [10].

Theorem 1** (Friedlander et al., 1999).**

Consider the linear system (1) with $A\in\mathbb{R}^{n\times n}$ is SPD and $b\in\mathbb{R}^{n}$, where $x_{*}=A^{-1}b$ is the exact solution. Consider the gradient method (3) being used to solve (1) and the steplength $\alpha_{k}$ given by (10). Then the sequence $\{x_{k}\}$ converges to $x_{*}$ starting from any point $x_{0}$.

For a proof of the above theorem, see [10]. Incidentally, several potential algorithms were provided therein, including the first cyclic gradient method under the name of cyclic steepest descent (CSD) as suggested in [4], which can be summarized as follows

[TABLE]

Notice that if we choose $\rho(k)=0$ and $\tau(k)=\bar{k}+1,\ \dots,\ k-1,\ k$, then (10) becomes CSD method, which satisfies the Theorem 1. On the other hand, Dai [4] proposed a variant called cyclic Barzilai-Borwein (CBB) method. They suggested that

[TABLE]

Similarly, if we choose $\rho(k)=0$ and $\tau(k)=\bar{k},\ \bar{k}+1,\dots,\ k-1$, then (10) becomes CBB method.

Although these methods greatly speed up the convergence, their motivation is too straightforward to further accelerate the iterations, which relies on the nonmonotone property to search the whole space without sink into any lower subspace spanned by eigenvectors [9]. This allows to reduce the gradient components more or less in the same asymptotic rate [6].

The recent literature showed that Yuan steplength may lead to efficient algorithms [20, 6]. All methods therein have two-dimensional finite termination property, i.e., if (8) is applied to a linear system in two-dimensional space, then the algorithm will terminate in at most 3 iterations. In general, such property seems not attractive in practice. However, experiments showed that they perform well in higher dimensions and are competitive with BB methods for large-scale problems [6].

Inspired by the Yuan steplength, we suggest a simple way of modifying steepest descent model to a cyclic gradient method. Consider a steplength of the form

[TABLE]

Here we modify the order of SD and Y compared to the original YB formula, which is useful for the development of the new algorithm. Apart from this change, (15) is indeed the second algorithm propose by the pioneering work of Yuan [20]. It keeps the two-dimensional finite termination property that performs as well as BB for large-scale problems and better for small-scale problems. We could introduce simply the cyclic behavior based on (15) of the form

[TABLE]

Besides, we find that De Asmundis et al. [8] gives an interesting view about the iterations of SD method, where the technique of alignment was proposed therein to force the gradients into one-dimensional subspace and avoid the zigzag pattern. Notice that the inverse of constant Rayleigh quotient such as SD and BB steplengths has the similar property. Thus, constant SD with retards can also give rise to the alignment behavior and keep the nonmonotone benefit. To achieve this goal, we need to impose a repeat time to the zigzag process. Meanwhile, we want to keep the process based on Yuan steplength in the first several iterations. These motivations lead to a new method of the form

[TABLE]

where $l\geq 1$ and $m\geq 1$. Such formula seems complicated, but indeed easy to understand. There are three components consisting in (17): the first SD and Y are used to insure the finite termination property; the parameter $l$ acting on the second part of SD is used to keep several zigzag iterations; finally, the retard term $m$ induces alignment and provides nonmonotone behavior to leap from the lower subspace.

3 Convergence Analysis

By the invariance property under any orthogonal transformation, we can assume without loss of generality that

[TABLE]

where

[TABLE]

We follow the convergence framework established by Dai [4] and adapt it to our method. Let

[TABLE]

where $g_{i,k}$ is the $i$th component of $g_{k}$. A preliminary property is defined as follows.

Definition 2** (Property A).**

Suppose that matrix $A$ has the form (18) with condition (19) holds. If $\exists\xi\in\mathbb{N}$, $\exists M_{1},M_{2}>0$, such that $\forall\mu\in\{1,\dots,n-1\}$, $\forall\epsilon>0$, $\forall j\in\{0,\dots,\min\{k,\xi\}-1\}$,

• •

$\lambda_{1}\leq\alpha_{k}^{-1}\leq M_{1}$;

• •

if $G(k-j,\mu)\leq\epsilon$ and $g_{\mu+1,k-j}^{2}\geq M_{2}\epsilon$, then $\alpha_{k}^{-1}\geq\frac{2}{3}\lambda_{\mu+1}$,

then the steplength $\alpha_{k}$ has Property A.

The convergence framework of Dai can be deduced from Property A, stated as follows.

Theorem 3** (Dai, 2003).**

Consider the linear system (1) with $A\in\mathbb{R}^{n\times n}$ of the form (18) and $b\in\mathbb{R}^{n}$. Consider the gradient method (3) being used to solve (1). If the steplength $\alpha_{k}$ has Property A, then the sequence $\{\left\lVert g_{k}\right\rVert\}$ converges to 0 R-linearly for any starting point $x_{0}$.

For a proof of the above theorem, see [4]. Many gradient methods have Property A as mentioned in [4], e.g., the gradient method with retards (10). Inspired by the demonstration therein, we now develop a convergence result for the CY method.

Theorem 4**.**

Consider the linear system (1) with $A\in\mathbb{R}^{n\times n}$ of the form (18) and $b\in\mathbb{R}^{n}$. Consider the gradient method (3) with steplength (17) being used to solve (1). Then the steplength $\alpha_{k}^{\text{CY}}$ has Property A.

Proof.

Note that (17) has three alternate steplengths, whereas the SD updating process and the constant process using the last SD steplength both follow the framework (10), which has been proven to have the Property A [4]. Therefore, we only investigate the Yuan steplength.

Recall that Yuan steplength has the following property

[TABLE]

which is given in [20]. Hence,

[TABLE]

Then the first condition of Property A holds by setting $M_{1}=2\lambda_{n}$. For the second one, let $M_{2}=2$ and $\xi=1$, which yields $j=0$. Suppose that

[TABLE]

for all $\mu\in\{1,\dots,n-1\}$, and $\epsilon>0$. Hence, the inverse of Yuan steplength becomes

[TABLE]

Hence, the second condition of Property A is satisfied, which completes the proof. ∎

4 Numerical Results

We first address the issue of parallel implementation. The dot product is engaged in the computation of steplength, which is the major obstacle of parallelization. Here we have two strategies to realize this goal. Let $A_{i}$ be the band matrix stored in the $i$th processor. The first one (Gather Algorithm, GA) is to gather the vector $q_{i}=A_{i}*g$ and execute dot product with global vectors, shown as follows

  $\texttt{Allgatherv}(q,\ q_{i})$

  $\alpha=\texttt{Dot}(g,\ g)\ /\ \texttt{Dot}(g,\ q)$

Then, the second one (Reduce Algorithm, RA) consists in computing the dot product locally, shown as follows

  $c_{i}=\texttt{Dot}(g_{i},\ q_{i})$

  $\texttt{Allreduce}(c,\ c_{i},\ \texttt{SUM})$

  $\alpha=\texttt{Dot}(g,\ g)\ /\ c$

Besides, we can see that global gradient vector is used in each iteration that we must proceed another Allgatherv function to communicate with other processor.

Let $p$ be the number of processors. The two experiments are proceeded by Alinea [12] (see also, e.g., [17, 18, 16]) and JACK [13, 15] (see also, e.g., [14]) and results are illustrated in Figures 1 and 2. We can see that generally the results are not good because the first one imposes so much computation and communication load, while the second one causes indeed the problem of loss of precision. These problems exist in all projection methods and by now we have not yet managed to find a solution.

The second experiments are proceeded by Matlab R2017b with a large-scale problem provided by The SuiteSparse Matrix Collection [7], with $n=50000$ and $349968$ non-zero values. All parameters are chosen under a training problem given in [2], such that $l=4,m=3$ for CY, $m=3$ for CSD, and $m=4$ for CBB. Average results are shown in Table 1.

We use bold numbers indicating the most efficient algorithms under each residual threshold. Backslash represents a number of iterations more than $10000$. From Table 1, we see that the CY method performs better than other methods except CG. Although we do not yet beat CG in high precision, our method is still competitive in most cases. Specifically, we can see that CY is stable throughout the iterations and much better than CG when low accuracy is required.

5 Conclusions

In this paper, we have proposed a new gradient method and shown that it is very competitive with CG method and better than the others for large-scale problems. However, it is still lack of theoretical evidence supporting such results. It is also important to find a better parallelization strategy. Therefore it still remains to study the properties and high performance implementation of gradient methods.

Acknowledgment

This work was supported by the French national programme LEFE/INSU and the project ADOM (Méthodes de décomposition de domaine asynchrones) of the French National Research Agency (ANR).