On C-To-R-Based Iteration Methods for a Class of Complex Symmetric Weakly Nonlinear Equations

To avoid solving the complex systems, we first rewrite the complex-valued nonlinear system to real-valued form (C-to-R) equivalently. Then, based on separable property of the linear and the nonlinear terms, we present a C-to-R-based Picard iteration method and a nonlinear C-to-R-based splitting (NC-to-R) iteration method for solving a class of large sparse and complex symmetric weakly nonlinear equations. At each inner process iterative step of the new methods, one only needs to solve the real subsystems with the same symmetric positive and definite coefficient matrix. Therefore, the computational workloads and computational storage will be saved in actual implements. The conditions for guaranteeing the local convergence are studied in detail. The quasi-optimal parameters are also proposed for both the C-to-R-based Picard iteration method and the NC-to-R iteration method. Numerical experiments are performed to show the efficiency of the new methods.


Introduction
We consider the iterative solutions of nonlinear system of equations in the following form, where A = W + iT ∈ C n×n is a large, sparse, complex symmetric matrix, with W ∈ R n×n and T ∈ R n×n being the real parts and the imaginary parts of the coefficient matrix A, respectively.
Here, we assume that W and T are both symmetric positive and semidefinite (SPSD) and at least one of them being symmetric positive and definite (SPD). The right hand vector function φ : D ⊂ C n → C n is a continuously differential function defined on the open convex domain D in the n-dimensional C n . u ∈ C n is an unknown vector. When the linear term Au is strongly dominant over the nonlinear term By substituting u = x + iy, where x ∈ R n and y ∈ R n , we can rewrite the system of nonlinear Equations (1) as (Wx − Ty) + i(Tx + Wy) = R(u) + iI(u), where R(u) = real(φ(u)) ∈ R n and I(u) = imag(φ(u)) ∈ R n are the real parts and imaginary parts of φ(u), respectively. Here, we reformulate the complex nonlinear system (1) as a block two-by-two real form in the following, When R(u) and I(u) are both constant vectors, the system (2) reduces to a linear system with structural block two-by-two coefficient matrix. To solve the linear system (2) efficiently, many methods have been proposed. Particularly, when the block matrices W and T are both SPSD and at least one of them is SPD, some classical iterative methods can be found in existed references, for example, the block preconditioned methods [9,10], the additive block diagonal preconditioned method [11]. There are also sequences of preconditioners based on the modified Hermitian and skew-Hermitian splitting (MHSS), such as, the preconditioned MHSS (PMHSS) iteration method [12,13]. Some other methods such as the preconditioned GSOR (PGSOR) iteration method [14], the complex-value to real-value (C-to-R) preconditioner [15], and so on, also attract a lot of researchers' interest. For more efficient methods, we refer to the works in [16][17][18].
When Au is strong dominant than the nonlinear term φ(u), and R(u) and I(u) are dependent on the variable vector u, then we say the system (2) is weakly nonlinear. The most classic and important solvers for the system of nonlinear Equations (1) is the Newton method [3,6,19], which can be described as u (k+1) = u (k) − F (u (k) ) −1 F(u (k) ), k = 0, 1, 2, · · · .
It can be seen from the Newton method that the dominant task in implementations is to solve the following equation at each iteration step, and to recompute the Jacobian matrix F (u (k) ) at every iteration step. When the Jacobian matrix F (u (k) ) is large and sparse, we usually use either the splitting relaxation form [6] or the Krylov subspace method form [5,20,21] to compute an approximation to update the vector s (k) . However, those methods are all heavy in both computational workload and computational storage.
To improve the efficiency of the Newton iteration method, Bai and Guo [4] use the Hermitian and skew-Hermitian splitting (HSS) method to solve approximately the Newton Equations (2), called the Newton-HSS method and then Guo and Duff [22] analyze the Kantorovich-type semilocal convergence. Many efficient methods based on the Hermitian and skew-Hermitian splitting (HSS) have come forth since then. For example, Bai and Yang [1] present the nonlinear HSS-like iteration method based on the Hermitian and skew-Hermitian (HS) splitting of the non-Hermitian coefficient matrix of the linear term Au. Some variants of the HSS-based methods for nonlinear equations can be found in references, e.g., the lopsided preconditioned modified HSS (LPMHSS) iteration method [23], the Newton-MHSS method [24], the accelerated Newton-GPSS iteration method [25], the preconditioned modified Newton-MHSS method [26], the modified Newton-SHSS method [27], the modified Newton-DPMHSS method [28], and so on. See [4,[29][30][31][32] for more details.
In this paper, we will concentrate on the efficient methods of the real value equivalent nonlinear system (2). By utilizing the C-to-R preconditioning technique proposed in [16], we will first construct a C-to-R-based Picard iteration method. This method is actually the inexact Picard method with the C-to-R iterative method as inner process iteration. Therefore, the convergence results based on Picard iteration method can be used directly. To further improve the efficiency, we then introduce a nonlinear C-to-R splitting iteration method for solving the weakly nonlinear equations. The local convergence for both of the two new methods are analyzed in detail. The way to determine the theoretical optimal parameters is also studied.
The organization of this paper is outlined as follows. In Section 2, we firstly give the C-to-R-based Picard iteration method. Then, we give the convergence results and the choice of the quasi-optimal parameters. In Section 3, we construct the nonlinear C-to-R-based splitting iteration method for solving the nonlinear system (2) and then we give a detailed theoretically analysis about the convergence properties. Theoretical optimal parameters are also proposed subsequently. Numerical experiments are given in Section 4 to illustrate the feasibility and effectiveness of the new methods. Finally, a brief conclusion and some remarks are drawn in Section 5 to end this work.
Throughout this paper, we use ρ(B) to denote the spectral radius of the matrix B. Denote u = x + iy ∈ C n with x ∈ R n and y ∈ R n being its real parts and imaginary parts, respectively.

The C-To-R-Based Picard Iteration Method
In this section, we will propose the C-to-R-based Picard iteration method. To begin with, according to [16], we review the C-to-R preconditioner for the block two-by-two coefficient matrix in Equation (2) as where α is a positive constant. We know that the implementing of the preconditioner B(α) at each iterative step needs one to solve the generalized residual linear equations Or equivalently, Therefore, we can summarize the implementation as the algorithm 1 for solving the above equations.
Algorithm 1 (The C-to-R iteration method). Let α be a given positive constant. Use the following steps to solve the generalized residual equation.
Step 2. solve (αW The preconditioner B(α) can be seen as a splitting matrix from the following matrix splitting, On the other hand, because Au is strong dominant over the nonlinear term φ(u), then the Picard iteration method can be used based on the separable property of the linear term and nonlinear term, i.e., Or equivalently, the iteration of the real value form can be described as where u (k) = x (k) + iy (k) ∈ C n with x (k) , y (k) ∈ R n . Therefore, to improve the efficiency, we can use the splitting iterative method as inner process iteration based on (6) to obtain the solution approximately at each Picard iteration process. We describe the C-to-R-based Picard iteration method in the Algorithm 2.
Algorithm 2 (The C-to-R-based Picard iteration method). Let φ ⊂ C n → C n be a continuously differentiable function and A = W + iT ∈ C n×n be a large, sparse complex symmetric matrix, where W ∈ R n×n and T ∈ R n×n are the real parts and the imaginary parts of A, respectively. Given an initial guess u (0) ∈ C n . For k = 0, 1, 2, · · · , until {u (k) } converges, compute the next iterate {u (k+1) } according to the following steps.
Ifz (k,l+1) does not meet the inner stopping criterion, return to Step 3.
From Algorithm 3, we find that the main workload is to solve the linear subsystem with the coefficient matrix being αW + T both in Step 4 and Step 5. Therefore, we can solve the corresponding system exactly by sparse Cholesky decomposition or inexactly by symmetric Krylov subspace methods (e.g., the conjugate gradient method).
Next, we denote then the local convergence results for the C-to-R-based Picard iteration method can be summarized in the following theorem.

Theorem 1. Denoted by
Then, for any initial guess u (0) ∈ C n and any sequence of positive integers l k , k = 0, 1, 2, . . ., if ν < 1 ln θ(α) , the iteration sequence generated by the C-to-R-based Picard method converges to the exact solution u and it holds lim or equivalently, the convergence rate of the C-to-R-based Picard method is R-rate, with the R-factor being ν.
Proof. The results can be obtained immediately from the results in [1].
Theorem 1 shows that the convergence rate of the C-to-R-based Picard iteration method depends on the quantities of θ(α) and ν. The weakly nonlinear property leads to the dominant of the quantities of θ(α). Therefore, we can obtain the quasi-optimal parameter by minimizing the R-factor θ(α) = L(α) . In other words, we need to find α such that the eigenvalues of the following matrix cluster around 1.
Therefore, we can use the results in [33,34] to obtain the quasi-optimal parameter α in the following theorem.

Remark 1.
Because there exists an extra term φ (u ), then we use the above optimal parameter as a quasi-optimal parameter. In actual implementation, we use this value as a suggestion. The true optimal parameter may vary dependent on the weakly nonlinear term.

The Nonlinear C-to-R-Based Splitting Iteration Method
In this section, we will further introduce the NC-to-R method for solving the block two-by-two nonlinear Equations (2) by making use of the nonlinear fixed-point equation The NC-to-R method can be described as the algorithm 4.
Algorithm 4 (The nonlinear C-to-R splitting iteration method). Let φ ⊂ C n → C n be a continuously differentiable function and A = W + iT ∈ C n×n be a large, sparse complex matrix, where W ∈ R n×n and T ∈ R n×n are the real parts and the imaginary parts of A, respectively. Given an initial guess u (0) ∈ C n .
The detailed implementing process of the NC-to-R method can be carried out in the Algorithm 5.
Next, we will focus on the convergence analysis for the NC-to-R method. By utilizing the Ostrowski Theorem (Theorem 10.1.3 in [6]), we can establish the local convergence theory for the NC-to-R method in the following theorem.
is defined in (2). If ρ( T(α; u )) < 1, then u ∈ D is a point of attraction of the NC-to-R iteration method.
After a few simple algebra computations, we have As φ : D ⊂ C n → C n is F-differentiable at a point u ∈ D, then the real parts and the imaginary parts are also F-differentiable at the point u . Therefore, Φ and Ψ are F-differentiable at the point u [35]. Therefore, by making use of the Ostrowski Theorem (Theorem 10.1.3 in [6]), we can conclude that if ρ( T(α; u )) < 1, then u is a point of attraction of the NC-to-R method.
At the end of this section, we will give a strategy to determine the optimal iterative parameter α by following the some strategy proposed in [36].
We use tr(.) to denote the trace of a matrix in the following. First, we do a partition of the matrix Φ (u ) as then the conjugate transpose of Φ (u ) can be expressed as where D 11 (.), D 12 (.), D 21 (.), and D 22 (.) have the same size as the matrix W and T.
The following theorem gives a strategy to choose the optimal parameter for the NC-to-R method. Then, the optimal parameter α opt for the NC-to-R method satisfies h (α opt ) = 0 and h (α opt ) > 0, Here, h (α) and h (α) are the first derivative and second derivative of h(α) with respect to the variable α, respectively.

Proof.
As it is known, if α is such a value such that B(α) is close to C, then B(α) −1 (C − B(α)) should be zero approximately. Therefore, C − B(α) should be approaching to zeros. Furthermore, the square norm C − B(α) + Φ (u ) 2 F could be reach the minimal with respect to α.

By direct algebra computations, it follows
Then, by using the notation in the theorem and after some simple algebra computations, we can obtain Therefore, by taking the first derivative of h(α) to be zero, it follows h (α opt ) = 0. Then the solution that satisfies the second derivative h (α) > 0 is the exact optimal parameter α opt that we need.

Numerical Experiments
In this section, we will testify the effectiveness of the C-to-R-based iteration methods by numerical experiments. All the tests are performed in MATLAB R2017a [version 9.2.0.538062] in double precision, on a personal computer with 2.40 GHz central processing unit (Intel(R) Core(TM) 2 Duo CPU), 4.00 GB memory and Windows 64-bit operating system. In our calculations, the stopping criterion for the proposed methods is the current relative residual satisfying r k 2 r 0 < 10 −6 , where r k is the residual at the k-th iteration with u (k) being the k-th approximate solution of Equation (1). We use the zero vector as the initial guess. To show the advantages of the new methods, we compare the C-to-R-based methods with the methods listed in Table 1, which shows the abbreviations and the corresponding full description. All the parameter choices in our experiments are listed in Table 2 and we classify the cases as follows. Table 1. Abbreviations and the corresponding description of the proposed tested methods.

N-G
Newton method using the GMRES method as inner iteration process N-H Newton method using the HSS method as inner iteration process HSS-N-G Newton method using the HSS preconditioned GMRES method as inner iteration process P-G Picard method using the GMRES method as inner iteration process P-H Picard method using the HSS method as inner iteration process HSS-P-G Picard method using the HSS preconditioned GMRES method as inner iteration process N-HSS nonlinear HSS-like method P-C Picard method using the C-to-R ietrative method as inner iteration process CP-P-G Picard method using the C-to-R preconditioned GMRES method as inner iteration process N-C nonlinear C-to-R-based iteration method In addition, we set the stopping criterion for the inner iteration process of all the methods to be where l k is the inner iteration steps number, η k is the prescribed inner tolerance. Here, we fix the η k simply by η = 0.1 for all k.
By applying the centered finite element different scheme with the space step size h = 1 N+1 and the implicit scheme with the temporal step size t = h, we can obtain the following nonlinear equations Here n = N 2 and ⊗ denoted the Kronecker product symbol. I n and I N are identity matrices of size n × n and N × N, respectively.  Our numerical results are presented for problem sizes N = 16, 32, 64, 128, 256, and 512 (i.e., n = 16 2 , 32 2 , 64 2 , 128 2 , 256 2 , and 512 2 ). First, we search the optimal parameter that minimizes the inner iteration count number by varying the parameter from 0.1 to 1 by step size 0.1. If the iteration number decreases when α decreases, then we will expand the searching area with extra [0.01, 0.1] by step size 0.01, or further [0.001, 0.01] by step size 0.001. Therefore, in Table 3, we give the experimental optimal parameters for all the proposed methods with respect to different cases and mesh grids.  The numerical results along with the it_out (i.e., the outer iteration counts), IT (i.e., the total inner iteration counts running through the corresponding method), and CPU (i.e., the elapsed cpu time in seconds) are shown in Tables 4-9. If the total inner iteration count number exceeds 500, or the elapsed CPU time is over 500 in second, or the computational storage is out of memory (especially the Newton-based methods), then we will denote the corresponding results as "-" in the tables.
From Tables 4-9, we find that the outer iteration counts of the Newton-based methods are less than the outer iteration counts of the Picard-based methods. However, the Newton-based methods occupy more CPU time than the Picard-based methods.   Besides, we find that the HSS-based iteration methods are more efficient than the Krylov subspace-based methods in CPU time. Further, the N-HSS method can solve the proposed problem efficiently with the optimal experimental parameters shown in Table 3, keeping the steady iteration count numbers while the mesh grid increases.
However, we also find that as the mesh grid increases, the CPU time of the N-HSS method increases rapidly.
The good news is that the C-to-R-based iterative methods (e.g., P-C and N-C) need the least iteration counts and CPU time. In particular, the NC-to-R method not only keeps a steady iteration count number with respect to different mesh grids, but also increases very slowly in CPU time as the mesh grid increases. Therefore, we can draw a conclusion that the C-to-R-based iterative methods are the first choice and best choice among all the proposed methods for solving this class of complex symmetric weakly nonlinear equations.

Concluding Remarks
In this paper, we focus on the numerical methods for solving a class of weakly nonlinear complex symmetric equations. First, we rewrite the original system as a real-valued form. Then, we propose a C-to-R-based Picard iteration method. This method is actually an inexact Picard iteration method. Therefore, the local convergence can be obtained by making use of some existed results. To further improve the efficiency, we construct a nonlinear C-to-R-based splitting iteration method. The convergence results and the theoretically optimal parameters are analyzed in detail. To illustrate the feasibility and the efficiency of the new methods, we perform some numerical experiments to compare with some classical methods. The numerical results show that our new methods are the most efficient method among all the proposed methods.