bugGNU Octave - Bugs: bug #44491, svd is slow and uses too much...

 
 

bug #44491: svd is slow and uses too much memory

Submitter:  Max Görner <maxg>
Submitted:  Mon 09 Mar 2015 07:12:06 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Performance
Status:  Need Info Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 3.8.2 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 18 May 2015 09:28:56 PM UTC, comment #8: 

I think we can close it, unless someone else still finds problems

Nir Krakauer <nir_krakauer>
Mon 18 May 2015 08:52:40 PM UTC, comment #7: 

@Nir: Is there any change to core Octave that needs to be done, or is this something exclusive to the statistics package?

Can we close this bug report?

Rik <rik5>
Group administrator
Tue 10 Mar 2015 07:05:02 AM UTC, comment #6: 

The patch works for me.

Max Görner <maxg>
Mon 09 Mar 2015 09:54:15 PM UTC, comment #5: 

I changed the script in the repository to use the economy SVD when n > p. This seems to solve the problem:

d = 0.01;
x = [-100:d:100]';
X = [x,x];
tic; s = princomp(X, "econ"); toc #took 9 s before fix, 0.002 s after

Nir Krakauer <nir_krakauer>
Mon 09 Mar 2015 03:35:30 PM UTC, comment #4: 

Max, to answer your question about copyright, since you've apparently read Matlab's source code, we cannot accept any code from you for this problem.

However, your suggestion to use an economy svd decomposition and/or implement Matlab's "0" argument for svd can be a good place to start.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Mon 09 Mar 2015 02:57:43 PM UTC, comment #3: 

From a previous IRC discussion, the problem is actually on svd being very very slow


2015-03-06 12:33:19        serviscopeMinor        well the super-super-super slow thing is confirmed to be in SVD. Not checked for results ...
2015-03-06 12:34:52        serviscopeMinor        I'll check for correctness now
2015-03-06 12:43:02        serviscopeMinor        s=100;
2015-03-06 12:43:02        serviscopeMinor        times=[];
2015-03-06 12:43:02        serviscopeMinor        r=0;
2015-03-06 12:43:02        serviscopeMinor        for i=1:100
2015-03-06 12:43:02        serviscopeMinor                n = floor(10000^(i/s));
2015-03-06 12:43:03        serviscopeMinor                data = rand(n, 64);
2015-03-06 12:43:05        serviscopeMinor                tic
2015-03-06 12:43:07        serviscopeMinor                a = svd(data');
2015-03-06 12:43:09        serviscopeMinor                t1= toc;
2015-03-06 12:43:11        serviscopeMinor                tic
2015-03-06 12:43:13        serviscopeMinor                [U, S, V]  = svd(data');
2015-03-06 12:43:15        serviscopeMinor                t2=toc;
2015-03-06 12:43:19        serviscopeMinor                r = r + sum(U(:)) + sum(V(:)) + sum(S(:));
2015-03-06 12:43:21        serviscopeMinor                errors = sum(abs(U*S*V' - data')(:));
2015-03-06 12:43:23        serviscopeMinor                times(end+1,:) = [n, errors, t1, t2];
2015-03-06 12:43:25        serviscopeMinor                i
2015-03-06 12:43:27        serviscopeMinor        end
2015-03-06 12:43:29        serviscopeMinor        hold off
2015-03-06 12:43:31        serviscopeMinor        clf
2015-03-06 12:43:33        serviscopeMinor        loglog(times(:,1), times(:,2), 'b-');
2015-03-06 12:43:35        serviscopeMinor        hold on
2015-03-06 12:43:37        serviscopeMinor        loglog(times(:,1), times(:,3), 'g-');
2015-03-06 12:43:39        serviscopeMinor        loglog(times(:,1), times(:,4), 'r-');
2015-03-06 12:43:41        serviscopeMinor        legend('errors', 'svd', '[u, v, s]= svd');
2015-03-06 12:43:43        serviscopeMinor        running that
2015-03-06 12:43:45        serviscopeMinor        no obvious errors, just massive slowness.
2015-03-06 12:47:24        serviscopeMinor        so there's certainly an awful bug in SVD.
2015-03-06 12:47:43        serviscopeMinor        in terms of speed
2015-03-06 12:49:12        serviscopeMinor        I'm just checking to see if all the SVDs are valid.
2015-03-06 12:49:21        serviscopeMinor        depending on the person's data, though the PCA is not unique.
2015-03-06 12:55:24        serviscopeMinor        carandraug: OK, the final script is this:
2015-03-06 12:55:31        serviscopeMinor        s=100;
2015-03-06 12:55:31        serviscopeMinor        times=[];
2015-03-06 12:55:31        serviscopeMinor        r=0;
2015-03-06 12:55:31        serviscopeMinor        for i=1:100
2015-03-06 12:55:31        serviscopeMinor                n = floor(10000^(i/s));
2015-03-06 12:55:32        serviscopeMinor                data = rand(n, 64);
2015-03-06 12:55:34        serviscopeMinor                tic
2015-03-06 12:55:36        serviscopeMinor                a = svd(data');
2015-03-06 12:55:38        serviscopeMinor                t1= toc;
2015-03-06 12:55:40        serviscopeMinor                tic
2015-03-06 12:55:42        serviscopeMinor                [U, S, V]  = svd(data');
2015-03-06 12:55:44        serviscopeMinor                t2=toc;
2015-03-06 12:55:46        serviscopeMinor                r = r + sum(U(:)) + sum(V(:)) + sum(S(:));
2015-03-06 12:55:50        serviscopeMinor                errors = sum(abs(U*S*V' - data')(:)) + sum(abs(U*U'-eye(size(U)))(:)) + sum(abs(V*V'-eye(size(V)))(:));
2015-03-06 12:55:52        serviscopeMinor                times(end+1,:) = [n, errors, t1, t2];
2015-03-06 12:55:54        serviscopeMinor                i
2015-03-06 12:55:56        serviscopeMinor        end
2015-03-06 12:55:58        serviscopeMinor        hold off
2015-03-06 12:56:00        serviscopeMinor        clf
2015-03-06 12:56:02        serviscopeMinor        loglog(times(:,1), times(:,2), 'b-');
2015-03-06 12:56:04        serviscopeMinor        hold on
2015-03-06 12:56:06        serviscopeMinor        loglog(times(:,1), times(:,3), 'g-');
2015-03-06 12:56:08        serviscopeMinor        loglog(times(:,1), times(:,4), 'r-');
2015-03-06 12:56:10        serviscopeMinor        legend('errors', 'svd', '[u, v, s]= svd');
2015-03-06 12:56:12        serviscopeMinor        so
2015-03-06 12:56:14        serviscopeMinor        apparently all the SVDs are valid
2015-03-06 12:56:16        serviscopeMinor        as in they reconstruct the original matrix and the matrices are all orthogonal.
2015-03-06 12:56:20        serviscopeMinor        I believe that's sufficient for it to be correct.
2015-03-06 12:56:22        serviscopeMinor        but very, very, VERY slow.


Carnë Draug <carandraug>
Group Member
Mon 09 Mar 2015 01:13:13 PM UTC, comment #2: 

Thanks for your bug report. Reassigning this bug report to the statistics package and adding maintainer in CC.

Mike Miller <mtmiller>
Group Member
Mon 09 Mar 2015 07:24:21 AM UTC, comment #1: 

I forgot to add, that in the described kind the error occours in 32bit only. Under 64bit Octave tries to eat up all memory and was killed by STRG+C or reboot before the call to princomp returned.

Max Görner <maxg>
Mon 09 Mar 2015 07:12:06 AM UTC, original submission:  

I'm developing some code under Matlab and Octave. However, princomp in Octave miserably.

My minimum working example is:

d = 0.01;
x = [-100:d:100]';
X = [x,x];
princomp(X);
error: out of memory or dimension too large for Octave's index type
error: called from:
error:   /usr/share/octave/packages/statistics-1.2.3/princomp.m at line 70, column 18


'X' is not even 1MiB in size, but I get this error. Matlab performs this code.

Furthermore, Matlab is way faster. If we change 'd' to 'd = 0.05' and perform the same code again, Octave takes 0.2 seconds and Matlab 0.001 seconds. On other matrices the ratio of ~100 reproduces.

I think I traced that down to the fact that Matlab ...
a) uses a special case with [coeff,~] = eig(X'*X)
b) calls 'svd(X,0)' whereas Octave calls 'svd(X)'.

If I enter the '0' in Octave's princomp function, it works. It is still a lot slower though.

Are there any reasons that the '0' is omittet? Are there any reasons that the special case using the eig function is not used? I'd create a patch if these questions are answered. Especially I'd like to know whether doing that call to eig would violate any copyrights.

sincerely,
Max Görner

Max Görner <maxg>

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by jordigh (Updated the item)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by mtmiller
  • -email is unavailable- added by maxg (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-05-18 rik5 Open/ClosedOpen Closed
    2015-05-18 rik5 StatusNone Need Info
    2015-03-09 carandraug CategoryOctave Package Octave Function
        Releaseother 3.8.2
        Summarystatistics package: princomp is slow and uses too much memory svd is slow and uses too much memory
    2015-03-09 jordigh Summarystatistics package: princomp is inferior to Matlab\'s statistics package: princomp is slow and uses too much memory
    2015-03-09 mtmiller CategoryOctave Function Octave Package
        Release3.8.2 other
        SummaryOctave\'s princomp is inferior to Matlab\'s statistics package: princomp is inferior to Matlab's
        Carbon-Copy- Added asnelt

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code