Heavy Ball Restarted CMRH Methods for Linear Systems

: The restarted CMRH method (changing minimal residual method based on the Hessenberg process) using fewer operations and storage is an alternative method to the restarted generalized minimal residual method (GMRES) method for linear systems. However, the traditional restarted CMRH method, which completely ignores the history information in the previous cycles, presents a slow speed of convergence. In this paper, we propose a heavy ball restarted CMRH method to remedy the slow convergence by bringing the previous approximation into the current search subspace. Numerical examples illustrate the effectiveness of the heavy ball restarted CMRH method


Introduction
In this paper, we are concerned with the CMRH method (changing minimal residual method based on the Hessenberg process) introduced in [1,2] for the solution of n × n linear system Ax = b.Given an initial approximation x 0 , and letting the initial residual r 0 = b − Ax 0 , the CMRH method is a Krylov subspace method based on the m-dimensional Krylov subspace overcome the slow convergence, there are many other effectiveness restarting technologies [14][15][16][17] which are designed to improve the simplest version of the restarted GMRES methods.Nevertheless, the research on the restarted CMRH method is rather scarce.One of the reasons is because the basis vectors generated by the Hessenberg process are not orthonormal.This leads to the CMRH residual vector being not orthogonal to the subspace K m (A, r 0 ) or A × K m (A, r 0 ).Thus, the whole augmented subspace containing the smaller Krylov subspace with Ritz vectors or harmonic Ritz vectors is not still a Krylov subspace.(See the subspace (2.4) of [17] for more details on this augmented subspace.)This means that the restarting strategy by including eigenvectors associated with the few smallest eigenvalues into the Krylov subspace [15][16][17] is not available in the restarted CMRH methods.In this paper, inspired by the locally optimal conjugate gradient (LOCG) methods for the eigenvalue problem [18][19][20][21][22][23] and the locally optimal and heavy ball GMRES methods for linear systems [14], we propose a heavy ball restarted CMRH method to keep the benefit and remedy the lost convergence speed of the traditional restarted CMRH method.For traditional restarted CMRH method (i.e., CMRH(m)), each CMRH cycle builds a Krylov subspace for computing the approximate solution of the cycle, which is used as the initial guess for the next cycle.As soon as the approximation is computed, the built Krylov subspace is thrown away.Nevertheless, in the heavy ball restarted CMRH method proposed in this paper, for salvaging the loss of the previous search space, we take the previous approximation into the current search to bring sufficient history information of the previous Krylov subspace.
The rest of this paper is organized as follows.In Section 2, we briefly review the Hessenberg process and the restarted CMRH method, then introduce the algorithmic framework and implementation details of the heavy ball restarted CMRH method.Numerical examples are given in Section 3 to show the convergence behavior of the improved method.Finally, we give our concluding remarks in Section 4.
Notation.Throughout this paper, R n×m is the set of all n × m real matrices, R n = R n×1 and R = R 1 .I n is the n × n identity matrix, and e (n) j is its jth column.The superscript ".T " takes the transpose only of a matrix or vector.For a vector u and a matrix A, u(j) is u's jth entry, u(i : j) is the vector of components u(i), . . ., u(j), and A(i, j) is A's (i, j)th entry.Notations • 1 , • 2 and • ∞ are the 1-norm, 2-norm, and ∞-norm of a vector or matrix, respectively.

The Hessenberg Process with Pivoting
The CMRH method [1] for the linear systems Ax = b is based on the Hessenberg process with pivoting as given in Algorithm 1 to construct a basis of the Krylov subspace K m (A, r 0 ).Given an initial guess x 0 ∈ R n , the recursively computed basis matrix L m = [ 1 , . . ., m ] ∈ R n×m and the upper Hessenberg matrix H m ∈ R m×m by Algorithm 1 satisfy where m ) T and PL m is a unit lower trapezoidal matrix with In particular, the initial residual vector , where β 0 = r 0 (i 0 ) as defined in Line 2 of Algorithm 1.
Here, L † m+1 is the pseudo-inverse of L m+1 .In fact, any left inverse of L m+1 will work, but we use the pseudo-inverse here for simplicity.We can state equivalently that x m = x 0 + L m y m , where Like the quasi-minimal residual (QMR) method [24], the CMRH method is a quasi-residual minimization method.

The Heavy Ball Restarted CMRH (HBCMRH) Method
To alleviate a heavy demand on memory, the restarted CMRH method (or CMRH(m) for short) is implemented by fixing m and repeating the m-step CMRH method with the current initial guess , where the kth and (k − 1)th CMRH cycles are indexed by the superscript "(k)" and "(k − 1)", respectively.By limiting the number m in the Hessenberg process, although CMRH(m) successfully remedies possibly heavy memory burden and computational cost, at the same time it sometimes converges slowly.One of the reasons is that all the Krylov subspaces built in the previous cycles are completely ignored.Motivated by the locally optimal and heavy ball CMRES method [14] which is proposed by including the approximation before the last to bring in more information in the previous cycles, we naturally develop a heavy ball restarted CMRH method denoted by HBCMRH(m) to make up the loss of previous search spaces.For each cycle, HBCMRH(m) starts an approximation x Next, we run the similar procedure to annihilate the m + 1 components p(1), . . ., p(m + 1) of the vector A ˆ m+1 , and then obtain ˜ m+2 and H m+1 ∈ R (m+2)×(m+1) as Lines 14-20 in Algorithm 2. Let L m+2 = [L m+1 , ˜ m+2 ].We have the relationship At the same time, P L m+2 keeps the structure of unit lower trapezoid.
The new approximation solution of HBCMRH(m) is x Here, L † m+2 is the pseudo-inverse of L m+2 .For any z ∈ K m (A, r 0 ) + span{d} can be expressed by z = L m+1 y for some y ∈ R m+1 .It is followed by (2) Thus, the new HBCMRH(m) approximation solution x We summarize the HBCMRH(m) method mentioned in this subsection in Algorithm 2. A few remarks regarding Algorithm 2 are in order: 1.In Algorithm 2, we only simply consider the case d = 0 and u = 0 in Line 11 and 18, respectively.
In fact, in the case d = 0 and u = 0, by a simple modification of the above process, the new approximation x Determine ) and p(m + 1) ↔ p(j 0 ).

Numerical Examples
In this section, we present some numerical examples to illustrate the convergence behavior of the HBCMRH(m) method (i.e., Algorithm 2) with the initial vector x 0 = [0, . . ., 0] T and m = 30.In demonstrating the quality of computed approximations, we monitor the normalized residual norms against the number of cycles.All our experiments were performed on a Windows 10 (64 bit) PC-Intel(R) Core(TM) i7-6700 CPU 3.40 GHz, 16 GB of RAM using MATLAB version 8.5 (R2015a) with machine epsilon 10 −16 in double precision floating point arithmetic.
Example 1.We first consider a n × n dense matrix appearing in [1] and [25] with n = 100, a i = 1 + i × ε, and ε = 10 −2 .The right hand side b = rand(n, 1) where rand is a MATLAB built-in function.In order to be fair in comparing algorithms, we tested the HBCMRH(m) method and the CMRH(m + 1) method.Recalling the remark of Algorithm 2, we know the CMRH(m + 1) method requires the same matrix-vector multiplications as HBCMRH(m), and it presents a better convergence behavior than CMRH(m).The normalized residual norms against the number of cycles of these two methods are collected in Figure 1, which clearly shows that the HBCMRH(m) method converges much faster than CMRH(m + 1).In fact, as shown in Figure 1, to reach about 10 −8 in normalized residual norms on this example, the HBCMRH(m) method takes 34 cycles, while the CMRH(m + 1) method is seen to need much more cycles to get there.Example 2. In this example, we consider a sparse matrix A and right hand side b which are extracted from raefsky1 taken from the University of Florida sparse matrix collection [26].In such a problem, A is not symmetric with order n = 3242 and contains 293,409 nonzero entries.Similarly, we compare the normalized residual norms of the HBCMRH(m) method with the CMRH(m + 1) method in Figure 2. Similar comments to the ones we made at the end of Example 1 are valid here as well.Example 3. In this example, we use the symmetric Hankle matrix appearing in [8] with elements

Number of cycles
and n = 1000.Let the right rand side b = A × ones(n, 1), where ones is a MATLAB built-in function.We compare the HBCMRH(m) method with the heavy ball restarted GMRES method denoted by HBGMRES(m) and compute the associated normalized residual norms in Figure 3.As shown in Figure 3, in such a case, the HBCMRH(m) and HBGMRES(m) method are competitive in the number of cycles, but less computation cost and storage requirement are needed in the HBCMRH(m) method.

Conclusions
In this paper, we proposed a heavy ball restarted CMRH method (HBCMRH(m) for short; i.e., Algorithm 2) for a linear system Ax = b.Compared to the traditional restarted CMRH method (i.e., CMRH(m)), x (k−1) 0 is built in the search space of the HBCMRH(m) method to compute the next approximation x (k) m for alleviating the loss of previous Krylov subspace.Therefore, one more matrix-vector multiplication of A is required in the HBCMRH(m) method.However, as shown in our numerical examples, HBCMRH(m) presents a better convergence behavior than CMRH(m + 1).
We have focused on the heavy ball restarted CMRH method.In fact, we can easily give the locally optimal restarted CMRH method as the locally optimal restarted GMRES method.We omit the details.In addition, while the HBCMRH(m) method has been developed for real linear systems, the algorithm can be rewritten to work for complex linear systems.This is done by simply replacing all R by C and each matrix/vector transpose by complex conjugate transpose.
previous approximation solution x (k−1) m previous cycle approximation solution x (k−1) m and seeks the next approximation solution x

.
The actual implementation is as follows.Let the vector d = x Recall that we have L m , H m , and p by the Hessenberg process satisfying AL m = L m+1 H m , where PL m is unit lower trapezoidal with P T = [e (n) p(1) , . . ., e (n) p(n) ].Now we let d subtract multiples 1 , . . ., m to annihilate the m components p(1), . . ., p(m) of the vector d to obtain a new vector d new as in Lines 8-10 in Algorithm 2. Suppose d new = 0. We determine

1 −) 1 − H m y 2 . 2 .Algorithm 2 1 :
H m+1 y 2 where H m+1 ∈ R (m+1)×(m+1) is obtained by deleting the last row of H m+1 .Similarly, if d = 0, then y (k) m is the optimal argument of min y∈R m β 0 e (m+1In comparison with the CMRH(m) method, the HBCMRH(m) method only requires one extra matrix-vector multiplication with A, but it represents a significant improvement in the speed of convergence, as shown in our numerical examples.Heavy ball restarted CMRH method (HBCMRH(m)) Choose an initial guess x (1) 0 ∈ R n and an integer m ≥ 1. 2: for k = 1, 2, . . .until convergence do 3: call Algorithm 1 to get p, L m+1 , and H m .
The CMRH approximate solution x m is obtained as x 0 + z m , where z m solves min