[Octave-bug-tracker] [bug #55727] Feature request: add GEJSV as an additional svd_drivers

[count]

Follow-up Comment #1, bug #55727 (project octave):

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)
    _______________________________________________________

Additional Item Attachment:

File name: svd-add-gejsv-v1.patch         Size:34 KB
    <https://savannah.gnu.org/file/svd-add-gejsv-v1.patch?file_id=46299>

File name: test_gejsv_all_cases.m         Size:1 KB
    <https://savannah.gnu.org/file/test_gejsv_all_cases.m?file_id=46300>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?55727>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/