I am interested in the iterative solution (preferably Krylov-type solvers) of a problem Ax=b, with x, b ∈ ℂn x 1 and A ∈ ℂn x n. 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:
- simple algorithm (i.e. without the need for orthogonalization during each step, as in GMRES) which converges well with a diagonal preconditioner.
- in which the estimates of xn are guaranteed to "monotonically decrease" the residual rn = b - Axn. By this I mean that, for example, ||rn||2 or rnHArn 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, 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, 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 AHAx=AHb, where the superscript H is hermitian (complex conjugate). Since A is invertible, then AHA is necessarily hermitian. And 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 AHA will be equal to the condition number of 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 AHAx=AHb 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 ||rn||/||b|| for solving my problem using the Conjugate Gradient method (CGM) of AHA, 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 A of shape 6930 x 6930.
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?