[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Switch svd driver default to gesdd?
From: |
Rik |
Subject: |
Re: Switch svd driver default to gesdd? |
Date: |
Wed, 28 Dec 2016 22:23:47 -0800 |
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