bugGNU Octave - Bugs: bug #41453, norm is really slow for sparse...

 
 

bug #41453: norm is really slow for sparse matrices

Submitter:  Mattias Linde <linde>
Submitted:  Tue 04 Feb 2014 11:01:53 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sat 30 Jan 2016 10:40:56 PM UTC, comment #24: 

Thanks for the patches, Massimiliano.  When I try to apply the patch, it complains that it is invalid.

When I apply it manually, it still fails to compile.  I think the problem may be that it is in liboctave/numeric, while the previous use was inclusion from libinterp/dldfcn.  Did you have similar problems?

(file #36214)

Lachlan Andrew <lachlan>
Wed 23 Apr 2014 08:00:29 AM UTC, comment #23: 

I have submitted a new version of the patch that uses SVD just if ARPACK is enabled and distinguishes real from complex matrices, invoking the appropriate ARPACK routine.

(file #31238)

Massimiliano Fasi <mfasi>
Fri 04 Apr 2014 08:34:46 AM UTC, comment #22: 

@Jordi: Thanks, now it makes sense. The idea of preparing an header file doesn't seem so good anymore either. On the other hand, any caveat about adding a method to check whether a sparse matrix is complex, if as I suspect there isn't any method to do this?

@Carlo: ok, I agree with your idea. It should be possible to do it using some preprocessor directives, as there is already a HAVE_ARPACK condition that can be used (now I also see why).

Massimiliano Fasi <mfasi>
Thu 03 Apr 2014 03:13:45 PM UTC, comment #21: 

As Jordi noted, ARPACK is currently an optional requirement.

I think the correct approach here would be to
make norm use the SVDS based computation of the norm
if it is available and fall back to the current slow algorithm if it is not.

c.

Carlo de Falco <cdf>
Group Member
Thu 03 Apr 2014 03:01:06 PM UTC, comment #20: 


> preparing a header file for eigs_base.cc. I don't know why that hasn't been done before, and why also libinterp/dldfcn/__eigs__.cc includes the file I am using as a .cc.


You have to understand the arrangement of Octave and why these functions are in dldfcn. First, it's possible to compile Octave without ARPACK. Sometimes people do this. Secondly, "dldfcn" means "dynamically loaded functions", i.e. dynamically linked. The code in this directory gets turned into oct files instead of being compiled into liboctave or libinterp. The idea here is to only load the libraries into RAM with dlopen when the user requests them (this is part of the reason why the first call to eigs is slower than later calls). And the reason it doesn't have a header file is because it's just meant to be a dynamically loaded oct file, somewhat "external" to Octave.

It may make sense to turn ARPACK into a core library like LAPACK, but it seems weird to make ARPACK a hard requirement just for two Octave functions: eigs, and svds.

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Thu 03 Apr 2014 09:28:02 AM UTC, comment #19: 

Are you sure ARPACK does not discriminate?

http://www.caam.rice.edu/software/ARPACK/UG/node70.html#SECTION001010000000000000000

http://www.caam.rice.edu/software/ARPACK/UG/node90.html#SECTION001020000000000000000

http://www.caam.rice.edu/software/ARPACK/UG/node108.html#SECTION001030000000000000000

It's good that you have a plan in mind.
Would you care about sharing your plan with us?

If it is too complicated to write it here, you could
use a wiki page for writing the plan and then post the link here.

c.

Carlo de Falco <cdf>
Group Member
Thu 03 Apr 2014 08:58:29 AM UTC, comment #18: 

Thank you! Well, I have in mind what should be done, and in my opinion there are just a couple of modifications that would be necessary to have a good sparse norm:

  • implementing a - one-line probably - method to check wether a sparse matrix is complex, equivalent to the is_complex_type () method you have suggested. It can't be used in such environment because I'm dealing with SparseMatrix and SparseComplexMatrix objects, not with octave_value. I cannot use the latter because I'm working in liboctave.


  • preparing a header file for eigs_base.cc. I don't know why that hasn't been done before, and why also libinterp/dldfcn/__eigs__.cc includes the file I am using as a .cc. Anyway, I don't think that using an header file instead would mess up something.


I don't think it would make sense to discriminate between hermitian and non-hermitian matrices as ARPACK doesn't provide different functions for these two classes of matrices.

Massimiliano Fasi <mfasi>
Thu 03 Apr 2014 08:17:43 AM UTC, comment #17: 

Massimiliano,

That is an improvement in the algorithm, and is going
in the right direction.

I agree though that including a .cc file is a really really really dirty hack.

I understand you are trying to produce something quick
to impress us during the GSoC selection period, but proposing
a well thought plan which shows you understand what the real
problem is here and what would be the correct way to solve
it would serve this purpose much better.

So before rushing into coding, I suggest you take some time
to discuss:

1) What is the best algorithm for computing the 2-norm of a SparseMatrix

2) What facilities the Octave numerical libraries offer that could help in implementing that algorithm.

3) What changes, if any, would be required to the Octave code-base in order to make those facilities available to oct-norm.cc

after clearing up these points you can go back to coding.

P.S. to discriminate between real and complex matrices, have a look at the "is_complex_type ()" method. Do you think it would also make sense to discriminate Hermitian and non Hermitian matrices?

Carlo de Falco <cdf>
Group Member
Wed 02 Apr 2014 07:45:35 PM UTC, comment #16: 

I have submitted a new version of the patch. It uses a function defined in eigs_base.cc to compute just one eigenvalue and it's indeed much faster than before, even faster than norm (full (A)). The preformances are now comparable with the ones of svd(A, 1) that remains a little bit faster (some experiments soon).

There are still some problems that should be adressed, in particular:

  • the fact that I am including a .cc file.


  • the use of a complex functions even if the input matrix is real, that makes the computation rather slower that it could be for real matrices.


(file #31108)

Massimiliano Fasi <mfasi>
Tue 01 Apr 2014 02:31:42 PM UTC, comment #15: 

Well, I think you are right in both cases. I was looking at the patch just as an improvement of norm (), but I see now that's not the way I should. Moreover, if the benchmark is svds (A, 1) there's no use in making comparisons with the version I submitted. For a 8000 x 8000 matrix I compute 16000 eigenvalues, that requires for sure much more time than computing just the largest magnitude eigenvalue as svds (A, 1) does (it's indeed even faster than norm (full(A))). Also the remark about full matrices is absolutely right, but this issue should be fixed for free using ARPACK.

Massimiliano Fasi <mfasi>
Tue 01 Apr 2014 01:13:09 PM UTC, comment #14: 

the main point is not that you EIG computes all eigenvalues, but that it allocates FULL matrices, which is definitely not something you want to do in general.

c.

Carlo de Falco <cdf>
Group Member
Tue 01 Apr 2014 12:41:56 PM UTC, comment #13: 

Your matrices are non very sparse, can you try again with


A = sprandn(10000, 10000, 0.01); tic, norm (A); toc


as I had suggested before?

also, consider that your results should be compared to:


A = sprandn(8000, 8000, 0.1); tic, svds (A,1,'L'); toc
Elapsed time is 13.895 seconds.


I still see a factor 100 ...



Carlo de Falco <cdf>
Group Member
Tue 01 Apr 2014 12:20:47 PM UTC, comment #12: 

I hadn't got over 1000 x 1000, but to answer your question I have run some experiments, whose results are


>> A = sprandn(10, 10, 0.1); tic, norm(A); toc
Elapsed time is 0.000300884 seconds.
>> A = sprandn(100, 100, 0.1); tic, norm(A); toc
Elapsed time is 0.00775719 seconds.
>> A = sprandn(1000, 1000, 0.1); tic, norm(A); toc
Elapsed time is 4.24463 seconds.
>> A = sprandn(2000, 2000, 0.1); tic, norm(A); toc
Elapsed time is 33.6406 seconds.
>> A = sprandn(4000, 4000, 0.1); tic, norm(A); toc
Elapsed time is 220.658 seconds.
>> A = sprandn(8000, 8000, 0.1); tic, norm(A); toc
Elapsed time is 1587.18 seconds.


for the patched version of the norm and


>> A = sprandn(10, 10, 0.1); tic, norm(A), toc
Elapsed time is 0.0279908 seconds.
>> A = sprandn(100, 100, 0.1); tic, norm(A), toc
Elapsed time is 3.31705 seconds.
>> A = sprandn(1000, 1000, 0.1); tic, norm(A), toc
Elapsed time is 333.631 seconds.
>> A = sprandn(2000, 2000, 0.1); tic, norm(A); toc
Elapsed time is 1433.76 seconds.


for the old version. I haven't tried larger matrices with the old norm as it would have lasted hours.

Massimiliano Fasi <mfasi>
Tue 01 Apr 2014 10:47:10 AM UTC, comment #11: 

how large was the matrix you used for testing?

try


A = sprandn (1000, 1000, .01);


or


A = sprandn (10000, 10000, .01);



c.

Carlo de Falco <cdf>
Group Member
Tue 01 Apr 2014 08:46:28 AM UTC, comment #10: 

Thank you for your feedback. Yes, I know that EIG computes all the eigenvalues of the matrix and is hence not the best choice, as we need just the largest one. Of course, that could be easily done with eigs. Nevertheless, we cannot call eigs from the file I am modifying, and I agree that the best choice would be using an ARPACK call, as you suggested. Reading a bit of source code, I think I've spotted the function I should use, that is defined in liboctave/numeric/eigs-base.cc, but I'm facing some troubles, as there are many parameters to call and I can't really understand how it works. I am working on it, but I don't know wether I 'll be able to get there in a reasonable amount of time.

I have tested my implementation with the current implementation in Octave. In particular, I compared


norm(A)


with


norm(full(A))


getting that the latter is up to five time faster than the former. Of course, that's not good, but it should be better than the factor 1000 we had at the beginning.

Massimiliano Fasi <mfasi>
Tue 01 Apr 2014 08:30:45 AM UTC, comment #9: 

Massimiliano,

Good work digging into Octave's source code, but I'm afraid you
are still heading the wrong direction.

Do you understand the difference between the Octave functions "eig" and "eigs"?

This bug report was about norm being slow for sparse matrices,
did you try to benchmark your implementation against the current
implementation in Octave?

c.

Carlo de Falco <cdf>
Group Member
Tue 01 Apr 2014 08:11:52 AM UTC, comment #8: 

I have submitted a new patch that should mitigate the issue. I have modified the standarnd norm function and used the EIG class to compute the singular values.

It could still be improved, using a direct call to ARPACK as suggested, but I am still working on it.

(file #31093)

Massimiliano Fasi <mfasi>
Mon 31 Mar 2014 08:58:19 AM UTC, comment #7: 

Thank you for the suggestion. I've explored a little the library and indeed it seems possible to avoid the use of svds (). The biggest problem now is concatenating matrices without using cat () or do_cat () from the core functions, but I hope I will find a solution.

Massimiliano Fasi <mfasi>
Sun 30 Mar 2014 06:12:59 PM UTC, comment #6: 

no it is not clean.

having to use a different
function to compute the norm of a sparse matrix
is not compatible and is no improvement at all
with having to call


svds (A, 1);


yes you can call an m-file function from c++
but that is not a clean solution either.

did you try to look at what svds does internally?
it will somehow rely on arpack in the end, and arpack
is a fortran library, not m-code.


do you think you would be able to implement this
using calls to arpack instead of svds.m?




Carlo de Falco <cdf>
Group Member
Sun 30 Mar 2014 05:41:53 PM UTC, comment #5: 

Ops, I submitter before finishing, sorry. Anyway, I was trying to say that I had introduced a new function because it seemed to me the cleanest way to do this. Do you think that "norm" could be fixed using the .m version of svds?

Thanks!

Massimiliano Fasi <mfasi>
Sun 30 Mar 2014 05:38:29 PM UTC, comment #4: 

I've tried to solve the issue introducing a feasible manner. I am not aware of a way to invoke an m function from a C++ file, that is the scenario I was facing. If such kind of trick can be done, then "norm" can be fixed easily. Otherwise that would require to rewrite svds in C++.

Massimiliano Fasi <mfasi>
Sun 30 Mar 2014 05:02:42 PM UTC, comment #3: 

I can't see the use of this patch.
What is the point in introducing a new function?
Please fix "norm" instead.
c.

Carlo de Falco <cdf>
Group Member
Sun 30 Mar 2014 04:54:18 PM UTC, comment #2: 

I have submitted a patch that tries to solve the bug by introducing a new function spnorm () that, as the name says, is supposed to be a norm for sparse matrices. It solves the performance issue with the sparse 2-norm using Carlo de Falco's idea in comment #1.

I've written an m-file instead of a C++ function because svds is already an octave script.

(file #31077)

Massimiliano Fasi <mfasi>
Tue 04 Feb 2014 11:36:06 AM UTC, comment #1: 

I don't know how "norm (A, 2)" is currently computed for
a sparse matrix, but probably we should use "svds (A, 1, 'L')"
instead, at least with your test this works better:



for s = [50 100 150]
  for i = 1:5
    A = sprand (s, s, 0.1);

    tic();
    n1 = svds (A, 1);
    t1 = toc();

    tic();
    n2 = norm (full (A));
    t2 = toc();

    printf (["Norm diff: %g, sparse time: %g,", ...
             "full matrix time: %g. Speed x: %g\n"],
            abs (n1-n2) / n1, t1, t2, t1/t2);
  end
end


Carlo de Falco <cdf>
Group Member
Tue 04 Feb 2014 11:01:53 AM UTC, original submission:  

In many cases I've encountered, it is much faster to do
norm( full( A ) ) 
compared to just norm(A) when A is sparse.

Short test script and diary are attached.
To avoid uploading testmatrices, the script uses random matrices, but the problem is consistent and shows up every time.

Have tested this with octave 3.2.4, 3.6.0 and dev sources as of today, similar results.

The speed difference I see is about a factor of 1000 and sometimes more.

Mattias Linde <linde>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36214:  bug_41453.v2.cset added by lachlan (4KiB - application/octet-stream)
file #31238:  sparse_norm.patch added by mfasi (4KiB - text/x-diff)
file #31108:  sparse_norm.patch added by mfasi (2KiB - text/x-diff)
file #31093:  sparse_norm.patch added by mfasi (2KiB - text/x-diff)
file #31077:  spnorm.patch added by mfasi (5KiB - text/x-diff)
file #30452:  slownorm.m added by linde (314B - text/x-objcsrc)
file #30453:  slownorm.txt added by linde (1KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Updated the item)
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by mfasi (Updated the item)
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by linde (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-08-02 siko1056 StatusNone Patch Submitted
    2016-01-30 lachlan Attached File- Added bug_41453.v2.cset, #36214
    2014-04-23 mfasi Attached File- Added sparse_norm.patch, #31238
    2014-04-02 mfasi Attached File- Added sparse_norm.patch, #31108
    2014-04-01 mfasi Attached File- Added sparse_norm.patch, #31093
    2014-03-30 mfasi Attached File- Added spnorm.patch, #31077
    2014-02-04 linde Attached File- Added slownorm.m, #30452
        Attached File- Added slownorm.txt, #30453

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code