normest1

[Marco Caliari]

Dear all,

I wrote the function normest1, which estimates the 1-norm of a matrix. It 
is based on Algorithm 2.4 of "A block algorithm for matrix 1-norm 
estimation, with an application to 1-norm pseudospectra", by N. J. Higham 
and F. Tisseur. The function call is the same as Matlab's normest1. I 
checked the implementation by trying to reproduce the tables included in 
the paper (see my testnormest1.m). The results are quite satisfactory. 
Matlab's normest1 is usually faster (in term of iterations), but produces 
much less exact results. It seems to me that the exit criteria were 
relaxed in favor of speed of "convergence".

About the similar function normest (2-norm estimator) which I contributed 
to write some time ago, now I have to say that I do not like any more the 
idea to fix the seed of the random generator. In fact, none of the 
previous tests (run several times and check whether you get the exact 
norm) would be possible. Moreover, if interested in an estimate, it is 
reasonable to run a couple of times (or even more) and take the maximum. 
On the other hand, fixing the seed is always possible before calling the 
function.

Finally, normest1 accepts a function 'afun = @(x) A*x' as input.
I think that normest could do the same (even if Matlab's normest does 
not). The most useful application that comes to my mind is the estimate of 
the 2-norm of A^m, which could be computed without forming A^m, by a 
function which does, essentially,

y = x; for i = 1:m, y = A * y;, endfor

I apologize for not using Mercurial. Please let me know if normest1.m can 
be accepted in Octave as is or if I need to do something else.

Best regards,

Marco

normest1.m
Description: Text document

testnormest1.m
Description: Text document