3
$\begingroup$

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:

  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 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.

CGM is bad in this case Figure 1: Clearly, CGM of AHA is bad in this case.

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

enter image description here 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?

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Breno
    Sep 13 at 17:48
  • $\begingroup$ 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...) $\endgroup$
    – Breno
    Sep 13 at 18:21
  • 1
    $\begingroup$ Nitpick: in CG it's the A-norm of the error that decreases monotonically, not that of the residual. $\endgroup$ Sep 13 at 21:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.