Re: norm of sparse matrices

[Jordi Gutiérrez Hermoso]

2011/11/25 Marco Caliari <address@hidden>:
> On Fri, 25 Nov 2011, Jordi Gutiérrez Hermoso wrote:
>
>> On 25 November 2011 06:10, c. <address@hidden> wrote:
>>>
>>> On 25 Nov 2011, at 11:00, Marco Caliari wrote:
>>>
>>>> Dear users,
>>>>
>>>> anybody knows which algorithm is used to compute the Euclidean norm of a
>>>> sparse matrix? An iterative method with a random initial value (such as in
>>>> normest)? Because I sometimes see different results on the same matrix:
>>>>
>>>> octave:39> S=sprand(5,5,0.5);
>>>> octave:40> norm(S)
>>>> ans =  1.08208354239579e+00
>>>> octave:41> norm(S)
>>>> ans =  1.08208354232114e+00
>>>> octave:42> norm(S)
>>>> ans =  1.08208354232114e+00
>>>>
>>>> (not easy to reproduce this behaviour). I'm using 3.4.1.
>>>>
>>>> Cheers,
>>>>
>>>> Marco
>>>
>>> Looking at the comments in the source code at
>>>
>>> http://hg.savannah.gnu.org/hgweb/octave/file/c3c8f513cf1f/liboctave/oct-norm.cc#l317
>>>
>>> http://hg.savannah.gnu.org/hgweb/octave/file/c3c8f513cf1f/liboctave/oct-norm.cc#l408
>>>
>>> it refers to "Higham's method" which I guess is the one described here
>>> http://www.springerlink.com/content/vh62v0354536p236/
>>>
>>> but I'm not an expert on the topic so I cannot help further.
>>
>> It appears this is correct. I have added the reference to the source:
>>
>>   http://hg.savannah.gnu.org/hgweb/octave/rev/7a756a7e145b
>>
>> If you have difficulty accessing the paper, I have a copy that I may
>> lend to anyone who asks.
>>
>> Thanks,
>> - Jordi G. H.
>
> I gave a look at the paper and at Higham's code pnorm.m, but I don't see any
> random initial guess. On the other hand, if look at this run
>
> octave:1> format long e
> octave:2> S=sprand(5,5,0.5);
> octave:3> norm(S)
> ans =  1.25336637903553e+00
> octave:4> norm(S)
> ans =  1.25336637903553e+00
> octave:5> norm(S)
> ans =  1.25336638124303e+00
> octave:6> norm(S)
> ans =  1.25336638124303e+00
> octave:7> norm(S)
> ans =  1.25336638124303e+00
> octave:8> norm(S)
> ans =  1.25336638124303e+00
> octave:9> norm(S)
> ans =  1.25336638124303e+00
> octave:10> norm(S)
> ans =  1.25336638124303e+00
>
> you see that the norm changes at the third computation and then it remains
> constant. This is not a random behaviour. Can anybody reproduce it? I can
> *almost* always reproduce it.

Yes, I was just able to reproduce that, and you're right, there's no
randomisation in the algorithm. This may point to a deeper error,
perhaps UB. I'm investigating with the debugger now...

- Jordi G. H.