[SciPy-dev] implement reorthogonalization method in gmres

[Darcoux Christine]

Dear Pauli,

Thank you for your quick answer with a pointer to the litterature.

The result in the paper you cite assumes that the matrix is
well-conditioned, and I'm not sure I can count on that in a
general-purpose code which many people will use.

The books "Solving nonlinear equations with Newton's method" and
"Iterative methods for linear and nonlinear equations" (some parts are
availlable from google book) by C. T. Kelley gives a good presentation
on the importance of orthogonality. It also show a simple 3x3 example
where paying attention to orthogonality is very important.

Checking for orthogonality and reorthogonalizing the way there author
of these books do it costs almost nothing (see the link to the code in
my previous post). Matlab's "official" GMRES code uses Householder
reflections to build the basis and gets orthogonality automatically,
but at double the cost of Gram-Schmidt. The Sandia Trilinos code uses
classical Gram-Schmidt for good parallel performance and
reorthogonalizes at every iteration, so also costs double. The way
proposed by Kelley costs nothing because the reorthogonalization
happens rarely, so I do not see any reason not to do it. Of course,
this could be added as an option so users will decide if they want to
use it.

Best Regards,

Christine

2010/2/24, Pauli Virtanen <pav at iki.fi>:
> ke, 2010-02-24 kello 23:43 +0100, nicky van foreest kirjoitti:
>> I am not at all a specialist on this, but I recall that Numerical
>> Recipes states that the Gram-Schmidt process gives terrible results,
>> and should be avoided nearly always.
>
> There's a balance here: more accurate methods tend to require more work,
> and the orthogonalization is probably the most expensive part of GMRES.
> So one should use a method that is good enough.
>
>> Instead, singular value
>> decomposition should be used. To avoid the problem below, would is be
>> an idea to replace the GM algo by SVD?
>
> One could also QR factorize using Householder transforms, which should
> be more accurate than GS, and is probably less expensive than SVD; or
> reorthogonalize.
>
> But as far as I understand, this would not really improve the
> performance of GMRES: as shown in the article I linked to (and also in
> others), the quality of the orthogonalization is not expected to matter
> for GMRES -- modified Gram-Schmid is good enough. The point is that
> breakdown of orthogonality happens simultaneously with the solution
> reaching its maximal accuracy, so it will not essentially affect
> convergence or accuracy. I'd like to have a reference or a test case at
> hand stating otherwise before changing things here.
>
> As far as I understand, however, the situation with orthogonality is
> different if you want to solve eigenvalue problems with Arnoldi stuff.
> But that's in the domain of ARPACK, and I would trust without looking
> that the ARPACK people know what they are doing.
>
> Cheers,
> Pauli
>
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>