Re: mgorth usage

[Ed Meyer]

On Mon, Sep 17, 2012 at 12:34 PM, JuanPi <address@hidden> wrote:

On Mon, Sep 17, 2012 at 9:32 PM, c. <address@hidden> wrote:
>
> On 17 Sep 2012, at 21:28, c. wrote:
>
>>
>> On 17 Sep 2012, at 18:29, Jordi Gutiérrez Hermoso wrote:
>>
>>> Gram-Schmidt is numerically unstable, so there's no reason for Octave
>>> to have it. Instead, use orth, which uses the svd.
>>
>> Gram-Schmidt is numerically unstable, that's why mgorth uses _modified_ Gram-Schmidt which is more stable.
>> THERE IS reason to have mgorth as it is a required step in any algorithm based on the Arnoldi process [1],
>> like, for example, GMRES [2].
>>
>>> HTH,
>>> - Jordi G. H.
>> c.
>>
>> [1] http://en.wikipedia.org/wiki/Arnoldi_iteration
>> [2] http://en.wikipedia.org/wiki/Generalized_minimal_residual_method
>
> jordi,
>
> reading your message again I see I had initially erroneously interpreted it
> as suggesting that Octave does not need mgorth. I now understand this was not
> the case, sorry for the noise. Hope the references on "mgorth"'s algorithm are
> of use to JP anyway.
>
> c.

Thanks again, all very useful!

--
JuanPi Carbajal
-----
"The bad economist pursues a small present good, which will be
followed by a great evil to come, while the true economist pursues a
great good to come, at the risk of a small present evil." - Frédéric
Bastiat
-----
http://ailab.ifi.uzh.ch/carbajal/

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave

I put a patch to mgorth on the patch tracker to fix a few problems with the documentation.

BTW the Householder QR method (as used in lapack) is more stable than MGS, although
more expensive. Both result in the same QR decomposition if the columns are independent;
MGS is less accurate as the columns become dependent, or as the condition number
rises. What this means for mgorth is that if x is nearly dependent on the columns of v
the resulting y will not be very orthogonal to v.  Using lapack QR to do the same thing
results in a much better y, as demonstrated by:

a = rand(10,5);
[q r] = qr(a); 
v = q(:,1:5);
x = q(:,5) + 10*eps*q(:,6);  # x is almost dependent on v
[y,h] = mgorth(x,v);
v'*y
[q r] = qr([v x]);
y2 = q(:,6);
v'*y2


-- 
Ed Meyer