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 (usv' - a, "fro")
ans = 4.2907e-14
octave:4> svd_driver ("gesdd");
octave:5> [u, s, v] = svd (a);
octave:6> norm (usv' - 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)
|
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 (u1s1v1' - a, "fro")
ans = 9.7938e-12
octave:7> norm (u2s2v2' - 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.
|
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.
|
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.
|
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.
|
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.
|
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.
|
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
?
|
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.
|
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).
|
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.
|
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
|
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?
|
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
|