A Breakdown-Free Block COCG Method for Complex Symmetric Linear Systems with Multiple Right-Hand Sides

The block conjugate orthogonal conjugate gradient method (BCOCG) is recognized as a common method to solve complex symmetric linear systems with multiple right-hand sides. However, breakdown always occurs if the right-hand sides are rank deficient. In this paper, based on the orthogonality conditions, we present a breakdown-free BCOCG algorithm with new parameter matrices to handle rank deficiency. To improve the spectral properties of coefficient matrix A, a precondition version of the breakdown-free BCOCG is proposed in detail. We also give the relative algorithms for the block conjugate A-orthogonal conjugate residual method. Numerical results illustrate that when breakdown occurs, the breakdown-free algorithms yield faster convergence than the non-breakdown-free algorithms.


Introduction
Consider the following complex symmetric linear system with multiple right-hand sides: with A ∈ C n×n non-Hermitian but symmetric (i.e., A = A H , A = A T ), X, B ∈ C n×p , and 1 ≤ p n. Such systems arise from many practical and important applications, for example, electromagnetic scattering, quantum chemistry, the complex Helmholtz equation, and so on [1,2].
Due to simple calculations and less information required, block Krylov subspace methods are always designed to solve system (1) efficiently [3,4]. Recently, Tadano et al. presented the block conjugate orthogonal conjugate gradient (BCOCG) method [5], which can exploit the symmetry of A naturally. The BCOCG is also deemed a natural generalization of the conjugate orthogonal conjugate gradient (COCG) method [6][7][8] for solving systems (1). Besides COCG-type methods, the COCR method described in [2,7,8] can also exploit the symmetry of A when p = 1 in (1). In [1], Gu et al. introduced a block version of the COCR method (BCOCR) by employing the residual orthonormalization technique. However, rank deficiency is a common problem that can lead block Krylov subspace methods to breakdown. The main reason is that the block search direction vectors may be linearly dependent on the existing basis by the increasing of the iteration number [9,10]. Consequently, some useless information will affect the accuracy of the solution and the numerical stability. BCOCG and BCOCR may also encounter such a problem, especially when p of (1) is larger, and we will show this phenomenon in the example section. Hence, it is valuable to solve the rank deficiency problem, and finally to enhance the numerical stability of BCOCG and BCOCR. Motivated by [10], in this paper, we propose a breakdown-free block COCG algorithm (BFBCOCG) and a breakdown-free block COCR algorithm (BFBCOCR) that can efficiently solve the rank deficiency problem of BCOCG and BCOCR, respectively.
The convergence rate of the Krylov subspace method depends on the spectral properties of the coefficient matrix, for example, the eigenvalue or singular value distribution, field of values, condition of the eigensystem, and so on. Fortunately, preconditioning can improve the spectral properties [11]. In this paper, we present the preconditioned version of BFBCOCG and BFBCOCR in detail.
The rest of this paper is organized as follows. In Section 2, based on the orthogonality conditions, we propose the BFBCOCG and BFBCOCR algorithms and their preconditioned variants with their new parameter matrices. Some numerical examples are listed in Section 3 to show the efficiency of our new algorithms. Finally, some conclusions are given in Section 4.

The Breakdown-Free Variants of BCOCG and BCOCR
In this section, we present our main algorithms, i.e., the breakdown-free block COCG algorithm (BFBCOCG) and the breakdown-free block COCR algorithm (BFBCOCR) in detail. We first introduce the block COCG (BCOCG) and block COCR (BCOCR) methods briefly and then give the derivation of BFBCOCG and BFBCOCR with their orthogonality properties, respectively. In the end, the preconditioned variants of BFBCOCG and BFBCOCR are also proposed, denoted by BFPBCOCG and BFPBCOCR, respectively. We use an underscore "_" to distinguish the breakdown-free and the non-breakdown-free algorithms.

BCOCG and BCOCR
Let X 0 be an initial approximation, and let X k+1 ∈ X 0 + K k+1 (A, R 0 ) be the k + 1th approximate solution of system (1). Hence, the recurrence relation of BCOCG and BCOCR is as follows: [1] The difference between BCOCG and BCOCR is the different calculation formulas of matrices α k and β k in (2), which are derived by applying the following orthogonality conditions: Different choices of L lead to different algorithms:

BFBCOCG and BFBCOCR
If the block size p is large, then the vectors of the block search direction will inevitably be linearly dependent on the increasing of the iteration number for BCOCG, hence rank deficiency occurs. In the following, in order to overcome this problem, we consider applying the breakdown-free strategy to BCOCG and propose the breakdown-free block COCG algorithm (BFBCOCG). The rationale of BFBCOCG is extracting an orthogonal basis P k+1 from the current searching space by using the operation orth(·). Thus, compared with (2), the new recurrence relation becomes Therefore, again using the orthogonality condition (3), we can get the Lemma 1. Here, we denote (3), thus R T j R k = 0 and P T j R k = 0. Then P T j AP k = 0 can be obtained by the second orthogonality condition in (3).
Similarly, the following Theorem 1 is obtained to update the parameters α k and β k . (3), the value of parameters α k and β k in the recurrence relation (4) can be obtained by solving the following equations:

Theorem 1. Under the orthogonality condition
Proof. From Lemma 1 and (4), we have the following two equations: and So, solving (6), we can easily get the α k . Pre-multiplying α T k to (7), then from the third equation of (4), we have From Lemma 1, we have R T k R k+1 = 0, and by the first equation of (5), one has α T k P T k U k = R T k P k . Thus the above equation can be rewritten as which can be used to update matrix β k .
The following Algorithm 1 is the breakdown-free block COCG.

Algorithm 1 Breakdown-free block COCG (BFBCOCG)
1. Given the initial guess X 0 ∈ C n×p and a tolerance tol; Similar to BFBCOCG, we can also easily get BFBCOCR by using L = AK k+1 (A, R 0 ) in the orthogonality condition (3). The following Algorithm 2 is the breakdown-free block COCR.

Algorithm 2 Breakdown-free block COCR (BFBCOCR)
1. Given the initial guess X 0 ∈ C n×p and a tolerance tol; Compute

BFPBCOCG and BFPBCOCR
As we all know, if the coefficient matrix has poor spectral properties, then the Krylov subspace methods may not robust, while a preconditioning strategy can make it better [11]. The trick is preconditioning (1) with a symmetric positive matrix M, which approximates to the inverse of matrix A; we get the following equivalent system: Let M = LL T be the Cholesky decomposition of M. Then system (9) is equivalent to We add a tilde "˜" on the variables derived from the new system. Then applying the BFBCOCG method and its recurrence relations (4) to (10), we have the orthogonality conditions It is easy to see The orthogonality condition (11) become R k ⊥K k (MA, MR 0 ) and AP k ⊥K k (MA, MR 0 ).
Under the recurrence relation (12) and the orthogonality condition (13), we can get the following Lemma 2 and Theorem 2 to update the matrices α k and β k . Here, we omit the proof because it is like the proof of Lemma 1 and Theorem 1. We denote Z k = MR k , U k = AP k .

Lemma 2.
For all 0 ≤ j < k, R T j Z k = 0, P T j R k = 0, and P T j AP k = 0.

Remark 1.
Under the preconditioned strategy, the relations of the block residuals are changed from orthogonal for BFBCOCG to M-orthogonal for BFPBCOCG. Here, two vectors x and y are M-orthogonal, meaning x⊥ M y, i.e., y H Mx = 0.

Theorem 2.
Under the orthogonality condition (13), the value of parameters α k and β k in the recurrence relations (12) can be obtained by solving the following equations: The following Algorithm 3 is the breakdown-free preconditioned block COCG algorithm.

Algorithm 3 Breakdown-free preconditioned block COCG (BFPBCOCG)
1. Given the initial guess X 0 ∈ C n×p and a tolerance tol; Change the orthogonality conditions (11) to the following conditions: The breakdown-free preconditioned block COCR (BFPBCOCR) can be deduced with the similar derivation of BFPBCOCG. Algorithm 4 shows the code of BFPBCOCR.

Algorithm 4 Breakdown-free preconditioned block COCR (BFPBCOCR)
1. Given the initial guess X 0 ∈ C n×p and a tolerance tol; At the end of this section, we will give the complexity for six algorithms. They are the block COCG algorithm, block COCR algorithm, breakdown-free block COCG algorithm, breakdown-free block COCR algorithm, breakdown-free preconditioned block COCG algorithm, and breakdown-free preconditioned block COCR algorithm, which are denoted by BCOCG, BCOCR, BFBCOCG, BFBCOCR, BFPBCOCG, BFPBCOCR, respectively. The pseudocodes of BCOCG and BCOCR are from [1]. Table 1 shows the main costs per cycle of the six algorithms. Here, we denote as "block vector" the matrix with size n × ×p, "bdp" the dot product number of two block vectors, "bmv" the product number of a matrix with n × p and a block vector, "bsaxpy" the number of two block vectors summed with one of the block vectors being from multiplying a block vector to a p × p matrix, "LE" the number of solving linear equations with a p × p coefficient matrix, and "bSC" the storage capacity of block vectors.
From Table 1, we can see the last four algorithms need one more dot product of two block vectors than the original two algorithms, i.e., BCOCG and BCOCR. For the product number of a matrix with a block vector, BFBCOCR and BFPBCOCR are both twice BFBCOCG and BFPBCOCG, respectively. This may result in more time to spend in BFBCOCR and BFPBCOCR, especially for a problem with a dense matrix. We reflect on the phenomenon in Example 1.
From Tables 2 and 3, we can see that for these rank-deficiency problems, the breakdown-free algorithms perform better than the non-breakdown-free algorithms in CPU time and iteration number. For the first type of problem, although the iterations of BFBCOCR are fewer than BFBCOCG, the CPU time is nearly double because the matrices are all dense and the product number of matrix and block vectors for BFBCOCR is one more than BFBCOCG per iteration. For the first two matrices of the second type of problem, the difference between BFBCOCG and BFBCOCR is not obvious, the main reason being that the matrices are sparse. However, for the third matrix helmdate_N160, only BFBCOCR can solve the problem. This indicates that for the matrices from the discretization of the Helmholtz equation, BFBCOCR performs the most robust compared to the other five algorithms.

Conclusions
In this paper, we presented a breakdown-free block conjugate orthogonal conjugate gradient algorithm for complex symmetric linear systems with multiple right-hand sides. Based on the orthogonality conditions, we gave its two new parameter matrices. The preconditioned version is also proposed in detail. At the same time, we also present the breakdown-free version for the block conjugate A-orthogonal conjugate residual method with its preconditioned version. From the numerical examples, we realized that when the right-hand sides are rank deficiency, our four new algorithms perform better than other algorithms. Moreover, for Helmholtz equation problems, BFBCOCR and BFPBCOCR show more stable behavior than BFBCOCG and BFPCOCG; while for dense matrices problems, BFBCOCG and BFPBCOCG converge faster than BFBCOCR and BFPBCOCR. However, there is still a lack of theoretical analysis for the advantages of BFBCOCG and BFBCOCR, and even of BFPBCOCG and BFPBCOCR. All of these require further investigations.

Conflicts of Interest:
The authors declare no conflict of interest.