bugGNU Octave - Bugs: bug #55727, Feature request: add GEJSV as an...

 
 

bug #55727: Feature request: add GEJSV as an additional svd_drivers

Submitted by:  Eddy <count>
Submitted on:  Sat 16 Feb 2019 10:24:23 PM UTC  
 
Category:  Octave Function Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  Patch Submitted Assigned to:  None
Originator Name:  Open/Closed:  Open
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

Tue 19 Mar 2019 05:11:57 PM UTC, comment #4: 

Update to the patch of the GEJSV driver (replace previous one), see attachment bug55727_gejsv_v2.patch.

Fixed few bugs in previous patch: about size of singular vector under some condition, and about the test case of GEJSV. Also it is now follow the C++ style guide tightly.

Now the patch is been tested for all combination of: different matrix sizes (tall, square, fat) with different number of output arguments, different data types (single, double, complex single, complex double) and full SVD v.s. economy SVD. Tested against --enable-address-sanitizer-flags.

Personally, I have used it for a while, so I think it is mature enough now.

To be more clear about what is inside, the changeset is consist of:

Estimated
lines-of-code    What for
202     declare FORTRAN functions
136     C++ interface to FORTRAN utilizer functions
129     compute optimal work space size for GEJSV
132     C++ interface to GEJSV
88      svd C++ class
rest    Fsvd () and documents.

The test script is also uploaded for convenient. Just run test_gejsv_all_cases.m.

(file #46580, file #46581, file #46582)

Eddy <count>
Sun 24 Feb 2019 08:58:14 PM UTC, comment #3: 

As an additional optional driver this should not do any harm.

The patch adds about 800 lines of code.  Some comments are included.  On a first sight it looks good.  I did not test it yet.

Kai Torben Ohlhus <siko1056>
Project Member
Sun 24 Feb 2019 02:57:20 PM UTC, comment #2: 

Add more time test results for svd drivers GEJSV, GESDD and GESVD (see attachment).

Looks good generally, i.e. have almost the same time cost as svd_driver GESVD when asked for left and right singular vectors. Severe slow down happens when matrix size (max of rows and cols) is close to or higher than 10^6.

(file #46349)

Eddy <count>
Mon 18 Feb 2019 02:14:52 AM UTC, comment #1: 

See attachment for the patch (with comments and tests) and test m-file.

The lack of workspace size query in [D,S]GEJSV make this patch significantly larger and involves a bunch of logic to compute it...

However, at the end of the day, GEJSV is shown to be very promising: it passes all the extreme tests appear in bug #55564.

Notably this one:

N = 26;
A = compan (fliplr ([1, 1 ./ cumprod(1:N)]));

a = A(1,:);
c = [1, -(1-sumsq(a)), -sumsq(a(1:end-1))];
s0 = sort(1-[roots(c); zeros(length(a)-2,1)], 'descend') .^0.5;

svd_driver('gesvd');
s1 = svd (A);

svd_driver('gejsv');
s2 = svd (A);

disp('   formula      gesvd        gejsv');
disp([s0 s1 s2]);

   formula      gesvd        gejsv
   6.0890e+26   6.0890e+26   6.0890e+26
   1.0000e+00   3.7617e+09   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   1.0000e+00   1.0000e+00
   1.0000e+00   6.6233e-01   1.0000e+00
   6.6233e-01   2.6583e-10   6.6233e-01

And the time cost is only slightly higher than gesvd.

tocs = @(s) fprintf('t = %7.3f sec: %s\n', toc(), s);

a = rand (2000);
svd_driver ('gesvd');
tic; [u1,s1,v1] = svd (a); tocs(svd_driver ());
svd_driver ('gejsv');
tic; [u2,s2,v2] = svd (a); tocs(svd_driver ());

% output
t =  54.570 sec: gesvd
t =  56.850 sec: gejsv

Memory footprint is larger than GESDD, but totally acceptable.

lwork = extra workspace size. measured in number of double type number.
econ = economy-sized decomposition.
Tested by e.g. svd_driver('gejsv'); [u,s,v]=svd(rand(100,10), 'econ'), with modification to source code to print lwork.

m x n \ algo        SDD        SDDecon        JSV        JSVecon
4x4                268        268        4292        4292
100 x 100        30700        30700        10400        20600
10 x 100        3310        770        7370        4490
100 x 10        3310        770        7370        4490
1000 x 10        32110        770        36170        4490
10000 x 10        320110        770        324170        20010

The large memory cost even for very small problem is due to the routine ormqr, somehow it tells me that it needs space at least 64*(64+1)=4160 for optimal performance.

If we are affordable to accept GESVD as default, why not pay a little more to accept GEJSV.

cheers

(file #46299, file #46300)

Eddy <count>
Sat 16 Feb 2019 10:24:23 PM UTC, original submission:  

Currently Octave supports two svd_drivers: GESVD and GESDD.

However, in some extreme cases, both of then are not satisfactory in terms of relative accuracy.

Such extreme cases can occur in research. e.g. several examples in bug #55564, e.g. total least square method needs the singular vector corresponding to the smallest singular value. See ref1 below for more.

The Jacobi method is known to be more accurate than both QR (GESVD) and Divide&Conquer (GESDD) method, see [ref1]. LAPACK includes such implementations named GEJSV and GESVJ. Accuracy and speed performance is analyzed in [ref2].

The GESVJ, is seems an internal routine for GEJSV, probably not intended for general use.

I'm writing a patch to make GEJSV available in Octave, hope that works.

[ref1]: Jacobi’s Method is More Accurate than QR
[ref2]: The Singular Value Decomposition: Anatomy of Optimizing an Algorithm for Extreme Scale

LAPACK manual GEJSV
LAPACK manual GESVJ

Eddy <count>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46580:  bug55727_gejsv_v2.patch added by count (36KiB - text/x-patch)
file #46581:  test_gejsv_all_cases.m added by count (3KiB - text/x-objcsrc)
file #46582:  svd_std_self_test.m added by count (1KiB - text/x-objcsrc)
file #46349:  blas_speed_survey_svd.pdf added by count (52KiB - application/pdf - time cost of svd drivers.)
file #46299:  svd-add-gejsv-v1.patch added by count (34KiB - text/x-patch - a patch to add gejsv svd driver)
file #46300:  test_gejsv_all_cases.m added by count (2KiB - text/x-objcsrc - a patch to add gejsv svd driver)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by count (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-03-19 count Attached File- => Added bug55727_gejsv_v2.patch, #46580
        Attached File- => Added test_gejsv_all_cases.m, #46581
        Attached File- => Added svd_std_self_test.m, #46582
    2019-02-24 siko1056 Severity3 - Normal => 1 - Wish
        Priority5 - Normal => 3 - Low
        StatusNone => Patch Submitted
    2019-02-24 count Attached File- => Added blas_speed_survey_svd.pdf, #46349
    2019-02-18 count Attached File- => Added svd-add-gejsv-v1.patch, #46299
        Attached File- => Added test_gejsv_all_cases.m, #46300

    Back to the top


    Powered by Savane 3.5