bugGNU Octave - Bugs: bug #29487, Feature request: svd with 3 output...

 
 

bug #29487: Feature request: svd with 3 output variables should use dgesdd

Submitter:  None
Submitted:  Fri 09 Apr 2010 12:30:04 AM UTC
   
 
Category:  Libraries Severity:  1 - Wish
Priority:  5 - Normal Item Group:  Feature Request
Status:  Fixed Assigned to:  highegg
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.3.51
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 31 Aug 2023 07:50:04 PM UTC, comment #14: 

The gesdd is much faster than gesvd. I also checked the gesdd_bad_matrix.dat matrix attached - it now seems to do better. Wondering if it's time to replace the default svd_driver...?

octave:1> X=rand(1000);
octave:2> svd_driver()
ans = gesvd
octave:3> tic;[U S V]=svd(X,'econ');toc
Elapsed time is 4.84571 seconds.
octave:4> svd_driver('gesdd')
octave:5> tic;[U S V]=svd(X,'econ');toc
Elapsed time is 0.465545 seconds.
octave:6>
octave:7> load gesdd_bad_matrix.dat
octave:8> svd_driver('gesvd')
octave:9> [U S V]=svd(a);
octave:10> svd_driver('gesdd')
octave:11> [U2 S2 V2]=svd(a);
octave:12> norm(a-U*S*V')
ans = 1.4619e-14
octave:13> norm(a-U2*S2*V2')
ans = 9.2040e-15

Anonymous
Mon 03 May 2010 12:24:13 PM UTC, comment #13: 

Attached is an ill-conditioned matrix that confuses the gesdd driver.

octave:1> load gesdd_bad_matrix.dat
octave:2> [u, s, v] = svd (a);
octave:3> norm (u*s*v' - a, "fro")
ans =  4.2907e-14
octave:4> svd_driver ("gesdd");
octave:5> [u, s, v] = svd (a);
octave:6> norm (u*s*v' - a, "fro")
ans =  70.878

So yes, it should be used more carefully. Given that SVD is used with ill-conditioned matrices quite often, it's reasonable to stay with gesvd as the default.

(file #20403)

Jaroslav Hajek <highegg>
Mon 03 May 2010 11:57:57 AM UTC, comment #12: 

I added svd_driver to Octave. Now I'm getting this:

octave:1> a = randn (1000);
octave:2> svd (1);
octave:3> tic; [u1, s1, v1] = svd (a); toc
Elapsed time is 10.1431 seconds.
octave:4> svd_driver ("gesdd")
ans = gesvd
octave:5> tic; [u2, s2, v2] = svd (a); toc
Elapsed time is 1.70709 seconds.
octave:6> norm (u1*s1*v1' - a, "fro")
ans =  9.7938e-12
octave:7> norm (u2*s2*v2' - a, "fro")
ans =  3.8047e-12

So it's true that gesdd can be much faster. Both drivers seem to be comparably accurate on normal data, though the singular vectors differ significantly due to lots of close singular values:

octave:8> norm (u1 - u2, "fro")
ans =  46.562

The default is still gesvd. Add svd_driver ("gesdd") in your .octaverc to enjoy the speed-up. Beware of possible problems for highly ill-conditioned, though I failed to produce a matrix that demonstrates this.

Jaroslav Hajek <highegg>
Thu 29 Apr 2010 05:37:33 PM UTC, comment #11: 

OK, given that users can screw themselves either way (pseudo variable or overloaded function) if they ask for speed over accuracy, then go ahead an use the variable if that seems most convenient.

John W. Eaton <jwe>
Group administrator
Mon 12 Apr 2010 01:16:33 PM UTC, comment #10: 

I understand, but I think you shouldn't apply this thinking too generally.

We have lots of config pseudo-variables for user preferences, because that's the generally best and most convenient way to deal with user preferences. All that I'm saying is that the choice of the driver will be more a user preference rather than an algorithm feature, because probably all svd-based algorithms will work with both drivers and simply inherit their properties (speed vs. robustness).

Making it a config variable is just a logical consequence of this premise. If you disagree with the premise, please explain why.

Jaroslav Hajek <highegg>
Mon 12 Apr 2010 01:03:29 PM UTC, comment #9: 

I would prefer to avoid adding more global "variables" to Octave, if at all possible.  Since this seems to be a case where you can achieve the same result with an overloaded function, I don't see the justification for adding the "variable".

Yes, I understand the problem with nargout == 0.

John W. Eaton <jwe>
Group administrator
Mon 12 Apr 2010 11:06:09 AM UTC, comment #8: 


In
http://www.cs.berkeley.edu/~ejr/tmp/lug.pdf
 page 12,  one read about several  new algorithms of
LAPACK 3 all of them are interesting.

In 4.9 stability of svd methods will be discussed.

Anonymous
Mon 12 Apr 2010 10:45:40 AM UTC, comment #7: 

I can't really agree, John, I think you're overgeneralizing.
99.9% of svd-based functions can't tell which driver they need, because it depends on the input data, and to reliably discover ill-conditioned matrices you generally need to svd them first. I believe this more or less holds for all svd-based functions in Octave's library, as well as OctaveForge packages and the rest of the world. That's where it differs from global variables that affect technical details, such as sparse_auto_mutate, which really is bad (and I'd vote for removing it, but that's a different matter).

So, if we provide an extra function, we won't be changing anything else to stay on the safe side, so users won't see any difference in svd-based functions, although they will be able to use svddc in their own scripts.

As you said yourself, most users (incl. me) will want to "just try" the new driver to possibly get the results faster. For their own scripts, both solutions are OK, but when they use other svd-based functions (such as the popular leasqr), svd_method is far more convenient than having to create @double/svd (and possibly also @single/svd) overloads.
Not speaking about the fact that it's a little tricky to write the overloads correctly (even you didn't do it quite right, watch out for nargout=0 case).

The solutions are interchangeable both ways, actually:
Just like you can override the library svd with svddc, you can also write svddc using unwind_protect + svd_method.
We can even distribute svddc with Octave, if you think it's worth the trouble.


Jaroslav Hajek <highegg>
Sat 10 Apr 2010 07:46:47 PM UTC, comment #6: 

If you make the change globally, then all functions that call svd will be affected, even those that need accuracy rather than speed.  So then we are in the situation of having to always write things like

  saved_svd_method = svd_method ();
  unwind_protect
    svd_method ("prefer_accuracy_over_speed");
    [u, s, v] = svd (m);
  unwind_protect_cleanup
    svd_method (saved_svd_method);
  end_unwind_protect

and I would need to remember to do this everywhere I called svd to avoid having some user's preference for speed over accuracy potentially screw up the calculation in the software I write and distribute (say, in some package).

I would rather not have to do that.  Exactly this kind of situtation is why I tried to remove all of the global "variables" that we used to have that could affect results of computations.

OTOH, I do understand that it is not convenient to have to change all the code that calls svd if you want to switch to using the other method globally, even just as an experiment to see what will happen.

What about using a different function name, say svddc (for svd divide-and-conquer, after the algorithm used in dgesdd).  Then if you really want to use this function globally in place of the default svd function, can't you do that by putting a @double/svd.m function in your path that does something like

  function varargout = svd (varargin)
    varargout = cell (nargout, 1);
    [varargout{:}] = svddc (varargin{:});
  endfunction

?

John W. Eaton <jwe>
Group administrator
Sat 10 Apr 2010 06:53:57 PM UTC, comment #5: 

I was thinking about a pseudovariable to control the driver. Do you see any problems with that? I agree that the default should be the current driver, i.e. GESVD, which should be more accurate for ill-conditioned problems.

Jaroslav Hajek <highegg>
Fri 09 Apr 2010 08:03:13 PM UTC, comment #4: 

How would you like to choose between the two methods?  Should it be a separate function?  Or an option to svd (we already have one optional argument, I think, so that complicates things slightly)?  I don't think you want to have a global change.  And I would definitely want the default to be the better/more accurate method (assuming one is clearly better than the other).

John W. Eaton <jwe>
Group administrator
Fri 09 Apr 2010 07:08:31 PM UTC, comment #3: 

Great. We can do the same, then. Actually I was not aware that GESDD is so much faster, what a shame. This could be an useful enhancement for me, I'll see what I can do.

Jaroslav Hajek <highegg>
Fri 09 Apr 2010 06:59:14 PM UTC, comment #2: 

It looks like the reverse issue came up in scipy which is claimed to use GESDD instead of GESVD and was producing incorrect results. Octave and Matlab were both cited as producing the correct results compared to scipy:

http://permalink.gmane.org/gmane.comp.python.scientific.devel/11002

It looks like scipy plans to make both available and let the user choose which engine to use:

http://projects.scipy.org/scipy/ticket/957

Judd Storrs <judd>
Fri 09 Apr 2010 08:19:43 AM UTC, comment #1: 

Is there any point in making this customizable? I.e., is GESDD generally superior to GESVD, or are there differences, for instance in accuracy?

Jaroslav Hajek <highegg>
Fri 09 Apr 2010 12:30:04 AM UTC, original submission:  

I recently checked the sage benchmark page
http://www.sagemath.org/tour-benchmarks.html
and found that matlab and octave are very slow generating
 [U s V]=svd(A)
the 1 output version of svd  is as fast as sage.

here my tests, on the same computer of course!:
octave:2> A=randn(1000);
octave:3> tic,[ U s V]=svd(A);toc
Elapsed time is 43.03 seconds.

sage:  m = random_matrix(RDF, 1000)
sage: %time U,s,Vh = m.SVD()
CPU times: user 6.38 s, sys: 0.28 s, total: 6.67 s
Wall time: 13.31 s
sage uses dgesdd.f

Ciau Erich

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #20403:  gesdd_bad_matrix.dat added by highegg (128KiB - chemical/x-mopac-input)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by judd (Posted a comment)
  • -email is unavailable- added by highegg (Posted a comment)
  •  

    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 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-05-03 highegg Attached File- Added gesdd_bad_matrix.dat, #20403
    2010-05-03 highegg Open/ClosedOpen Closed
    2010-05-03 highegg StatusNeed Info Fixed
    2010-04-29 jwe Severity3 - Normal 1 - Wish
    2010-04-09 highegg Assigned toNone highegg
    2010-04-09 highegg Item GroupNone Feature Request
        StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code