Re: Switch svd driver default to gesdd?

[Rik]

On 12/28/2016 01:42 PM, John W. Eaton wrote:
> On 12/28/2016 02:28 PM, Rik wrote:
>> 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.
>
> Thanks for checking.  I'm fine with switching the default if lapack
> generally gives better results now with gesdd.

I switched the default to gesdd on the development branch here
(http://hg.savannah.gnu.org/hgweb/octave/rev/750c8b4b7164).

Separately, I ran a few thousand random matrices through svd using both
drivers and checked accuracy with

[u,s,v] = svd (x);
norm (x - u*s*v', "fro")

Things have gotten better since 2010 because the gesdd driver had
consistently 50% smaller error than the gesvd driver.  The gesdd driver was
always faster and the performance multiple (Time (gesvd)   / Time (gesdd))
varied from 2.5X for small matrices to 8X for 4000x4000.

--Rik