Re: Switch svd driver default to gesdd?

[Rik]

On 12/28/2016 10:54 AM, Mike Miller wrote:
> On Wed, Dec 28, 2016 at 13:29:25 -0500, John W. Eaton wrote:
>> There was some discussion for this bug report:
>>
>>   https://savannah.gnu.org/bugs/?func=detailitem&item_id=29487
>>
>> Maybe some note of this should be added as a comment in the sources.
> Good digging. I wasn't even aware of the svd_driver function until now.
>
> I think it would be useful to add a note to the svd doc string to
> mention the speed / accuracy tradeoff and to explicitly point to
> svd_driver for more details (currently just listed as "see also").
>
> For comparison, scipy.linalg.svd takes the function name as an optional
> argument, so it is documented up front:
>
>     lapack_driver : {'gesdd', 'gesvd'}, optional
>         Whether to use the more efficient divide-and-conquer approach
>         (``'gesdd'``) or general rectangular approach (``'gesvd'``)
>         to compute the SVD. MATLAB and Octave use the ``'gesvd'`` approach.
>         Default is ``'gesdd'``.
>

I think Matlab must have switched over at some point to using gesdd.  The
bug report #29487 is from 2010 and Octave and Matlab were comparable in
speed, but slower than SAGE which had switched to gesdd.  Now, Matlab is
~40X faster than Octave as reported in bug #49940.  The issues raised in
the original bug report may also have been resolved.  An attachment to bug
#29487 is the file gesdd_bad_matrix.dat which produced bad results with
gesdd but acceptable results with gesvd.  Running that same file today
shows no significant difference so perhaps the LAPACK routine has been
improved.

--- Code ---
load gesdd_bad_matrix.dat
svd_driver ("gesvd");
[u, s, v] = svd (a);
norm (u*s*v' - a, "fro")

svd_driver ("gesdd");
[u, s, v] = svd (a);
norm (u*s*v' - a, "fro")
--- End Code ---

Results:
ans =    4.2078e-14
ans =    3.0277e-14

It still looks to me like we want to follow Matlab and use gesdd as the
default, but also add more documentation to svd and svd_driver.

--Rik