Re: Improved normest

[Marco Caliari]

Dear Jaroslav,

what about the following solution for the initial vector?

1) save the state of the random number generator with

v = rand("state");

2) set the state to something fixed (eventually depending on the matrix, 
such as

rand("state",trace(A))

3) "randomly" initialize the vector

x = rand(size(A,1),1);

4) perform the powers method and finally, before endfunction, restore the 
state of the random number generator with

rand("state",v).

Cheers,

Marco

On Thu, 26 Nov 2009, Jaroslav Hajek wrote:

> On Thu, Nov 26, 2009 at 2:09 PM, Marco Caliari <address@hidden>wrote:
>
> > On Thu, 26 Nov 2009, Jaroslav Hajek wrote:
> >
> >  On Thu, Nov 26, 2009 at 12:19 PM, Marco Caliari <address@hidden
> > > > wrote:
> > >
> > >  Dear maintainers,
> > > > I have seen changeset
> > > > http://hg.savannah.gnu.org/hgweb/octave/rev/87595d714005 about normest,
> > > > by
> > > > Jaroslav. As already discussed, the choice
> > > >
> > > > x = norm (A, "columns").'
> > > >
> > > > for the initial vector leads to wrong values (second largest eigenvalue)
> > > > for some "important" matrices, like
> > > >
> > > > toeplitz([-2,1,zeros(1,2)]), toeplitz([-2,1,zeros(1,4)]), etc
> > > >
> > > > This was the reason for the old random initial vector. I agree it is not
> > > > nice to have slightly different results at each run, but it is not nice
> > > > to
> > > > have wrong values, too.
> > > >
> > > > Best regards,
> > > >
> > > > Marco
> > > Already discussed? When? I couldn't find anything.
> >
> > Here, even if the choice for the initial vector was different
> >
> >
> > http://www-old.cae.wisc.edu/pipermail/octave-maintainers/2006-December/001422.html
> I see, thanks. I randomized the guess again. My problem with the previous
> code was that a priori forming A'*A is quite costly, especially if you want
> just low precision and expect just a dozen iterations:
>
> A = rand (1e3);
> tic; [e,c] = normest (A); toc
> tic; [e,c] = normest (A, 1e-1); toc
>
> now:
>
> Elapsed time is 0.0229969 seconds.
> Elapsed time is 0.0183249 seconds.
>
> previously:
>
> Elapsed time is 0.220653 seconds.
> Elapsed time is 0.210241 seconds.
>
>
>
> It's no problem to randomize it again, but I think a non-random behavior is
> > > preferable.
> >
> > I agree.
> Maybe there could be an option to specify the initial vector?
>
>
> >  Can the initial vector be improved somehow?
> > >
> >
> > I have no idea. It is usually written "take the initial vector you want:
> > roundoff errors will do the rest".
>
> Apparently reality is different :) But I have no idea either.
>
> --
> RNDr. Jaroslav Hajek
> computing expert & GNU Octave developer
> Aeronautical Research and Test Institute (VZLU)
> Prague, Czech Republic
> url: www.highegg.matfyz.cz