Inefficient handling of sparse matrices by the norm function

[David Bateman]

David Bateman wrote:
> Jason Riedy wrote:
>> And David Bateman writes:
>>> Yes, but still pathetic.. The issue is that the sum function is
>>> converting the matrix to full before doing the summation. I'm
>>> working on a patch
>> Oof.  I thought it was because I had a zillion things running.
>>
>> A quick patch might be to implement sum/spsum as a matrix-
>> vector product in the sparse case:
>>
>> octave:6> A = spdiag(rand(4E4, 1));
>> octave:7> tic; anorm = norm(A, 'inf'); toc
>> Elapsed time is 10.721107 seconds.
>> octave:8> tic; anorm2 = max(sum(abs(A),2)); toc
>> Elapsed time is 10.719147 seconds.
>> octave:9> tic; anorm3 = max(spsum(abs(A),2)); toc
>> Elapsed time is 10.686447 seconds.
>> octave:10> tic; anorm4 = max(abs(A)*ones(size(A,2),1)); toc
>> Elapsed time is 0.010323 seconds.
>>
>>
>> Jason
> 
> 
> Better to fix the src/data.cc{DATA_REDUCTION|NATIVE_REDUCTION} macros to
> accept sparse matrices.. Even better would be to add sum/cumsum, etc
> methods to the octave_value classes and use them rather than an ugly
> macro, but that can wait till after 3.0 :-)
> 
> I'm testing my patch now, as a update to 2.9.18+ caused a complete
> rebuild of Octave, so its a little slow..
> 
> D.


Well in fact it wasn't what I thought it was.. There was still some
unoptimized code for the sparse comparison, boolean and reduction
operators in Sparse-op-defs.h, but this affected not just "sum", but
things like "any" and "all", as well as the element by element
comparison and boolean operators. The attached patch addresses these
slowups. However, it doesn't modify the mixed sparse/matrix code that
still do dense indexing on the sparse matrices. However, as there is
already a dense matrix in the comparison or boolean operator, I don't
see that this will be a significant speed penalty.

I ran "make check" and all tests passed and so that appears to have not
introduced additional bugs. In addition I ran a number of additional
tests to confirm the correctness. Particularly nasty functions were
"prod" and "all" where the presence of the zeros had to be taken into
account (ie a column with a zero in it for all or prod results in zero,
where as pure sparse indexing would give a non-zero result).

The benchmark for the infinity norm of a sparse matrix is now.

A = spdiag(randn(1e4,1));
t0  = cputime();
B = norm (A, 'inf');
cputime () - t0

which gives the result of 0.048 seconds on my laptops where as it
previously took ~3 seconds. Note that the previous change in 2.9.17+ to
make the sparse to dense transformation optional and not the default,
means that norm(A,'inf') would returns, so I also forced the result to
be a scalar. There are probably other functions that take a sparse
matrix and return a scalar that have the same issue.

D.

-- 
David Bateman                                David.Bateman at motorola.com
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax)

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: patch1
Url: https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20071207/170c995d/attachment-0002.ksh 
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: patch1.changelog
Url: https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20071207/170c995d/attachment-0003.ksh