# Stable iterative solver for complex symmetric linear systems

I am interested in the iterative solution (preferably Krylov-type solvers) of a problem $$\boldsymbol{A}x=b$$, with $$x,b\in\mathbb{C}^{n\times1}$$ and $$\boldsymbol{A}\in\mathbb{C}^{n\times n}$$. $$\boldsymbol{A}$$ is symmetric, invertible, and its real and imaginary parts are indefinite. This comes from solving a harmonic, damped Finite Element problem.

Now, because the idea is to write a few iterations explicitly by hand, the challenge is that I would like to use a:

1. simple algorithm (i.e. without the need for orthogonalization during each step, as in GMRES) which converges well with a diagonal preconditioner.
2. in which the estimates of $$x_n$$ are guaranteed to "monotonically decrease" the residual $$r_n = b - \boldsymbol{A}x_n$$. By this I mean that, for example, $$||r_n||_2$$ or $$r_n^H\boldsymbol{A}r_n$$ are guaranteed to decrease monotonically per iteration. Sort of like the Conjugate Gradient method for hermitian matrices.

Are there any such methods?

I tried applying the Conjugate A-Orthogonal Conjugate Residual (COCR) method, by Sogabe and and Zhang, 2007[1], with simple diagonal preconditioning. I have a working implementation for coarse meshes, but for less coarse meshes, the algorithm fails to converge.

I then read about the Quasi-Minimal Residual COCR (QMRCOCR) method, by Gu et. al., 2014[2], which seems to be much more stable, still relatively simple to implement and the error seems to decrease monotonically. But this last remark is my interpretation of the results, not any rigorous convergence analysis. And a Symmetric Successive Over-Relaxation (SSOR) preconditioner is always used in the article, so it cannot be known whether the convergence will be as well-behaved if a simple diagonal preconditioner were used instead.

A possible solution that comes to mind is solving $$\boldsymbol{A}^H\boldsymbol{A}x=\boldsymbol{A}^Hb$$ where the superscript H denotes the hermitian adjoint. Since $$\boldsymbol{A}$$ is invertible, then $$\boldsymbol{A}^H\boldsymbol{A}$$ is necessarily hermitian (it is also positive-definite) Therefore, the simple Conjugate Gradient method could be used.

Now, everywhere I read about this mentions stability issues. I cannot find the reason why. The only concrete thing I have found is that the condition number of $$\boldsymbol{A}^H\boldsymbol{A}$$ will be equal to the condition number of $$\boldsymbol{A}$$ squared. But in my FE code, I always apply left- and right-diagonal preconditioning of A, thereby obtaining relatively good condition numbers.

So why exactly is solving $$\boldsymbol{A}^H\boldsymbol{A}x=\boldsymbol{A}^Hb$$ with an iterative solver bad?

After some experiments, the answer to the above question is: it just is.

Here are some images of the the curves $$||r_n||/||b||$$ for solving my problem using the Conjugate Gradient method (CGM) of $$\boldsymbol{A}^H\boldsymbol{A}$$, of the COCR method of the complex linear system and of the QMRCOCR method. I only used simple Jacobi preconditioning. It is a simple example for a small $$\boldsymbol{A}$$ of shape $$6930 \times 6930$$.

Figure 1: Clearly, CGM of $$\boldsymbol{A}^H\boldsymbol{A}$$ is bad in this case.

Figure 2: COCR and QMRCOCR both converge in 1132 and 1136 iterations, respectively. Clearly QMRCOCR is more stable.

Figure 3: The issue for me is that the residuals of neither COCR or QMRCOCR are guaranteed to decrease in the first iterations. For CGM it does decrease.

So is there a way to guarantee at least 2 first steps of a decreasing residual for COCR or QMRCOCR? I thought the first step of these types of methods was always similar to a gradient descent. Maybe for complex matrices, things work differently?

### References:

1. Sogabe, T.; Zhang, S. A COCR method for solving complex symmetric linear systems. J. Comput. Appl. Math. 2007, 199 (2), 297–303. DOI: 10.1016/j.cam.2005.07.032.
2. Gu, X.; Huang, T.; Li, L.; Li, H.; Sogabe, T.; Clemens, M. Quasi-Minimal Residual Variants of the COCG and COCR Methods for Complex Symmetric Linear Systems in Electromagnetic Simulations. IEEE Trans. Microwave Theory Techn. 2014, 62 (12), 2859–2867. DOI: 10.1109/TMTT.2014.2365472.
• After further study, it seems that the main issue is that diagonal preconditioning is simply not enough. Based on this Computational Science Stack Exchange answer by Daniel Shapero and based on this lecture by Wolfgang Bangerth, I cannot escape using a more complicated preconditioner while still having good convergence behaviors for an iterative solver. Sep 13 '21 at 17:48
• Then again, since A**<sup>H</sup>**A is guaranteed to be positive definite, then, at least for the initial steps of the secondary problem A**<sup>H</sup>**A*x*=**A**<sup>H</sup>*b*, the residual is guaranteed to decrease, right? (Edit: I cannot use HTML here...) Sep 13 '21 at 18:21
• Nitpick: in CG it's the A-norm of the error that decreases monotonically, not that of the residual. Sep 13 '21 at 21:35