bugGNU Octave - Bugs: bug #51880, sparse() ignores value of scalar...

 
 

bug #51880: sparse() ignores value of scalar column index input arg

Submitter:  Ray Zimmerman <rdzman>
Submitted:  Mon 28 Aug 2017 07:39:57 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Ray Zimmerman Open/Closed:  * Closed
Release:  * 4.2.1 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 06 Sep 2017 01:06:24 AM UTC, comment #8: 

Only way to know was to benchmark it.  Indeed, the resizing of the SparseRep was not a driver of execution time.  However, I did find that using std::fill_n over only the non-zero columns provided ~20% savings.

The benchmark code was


N = 1e3;

bm = zeros (N, 1);
for i = 1:N
  clear x;
  tic;
  x = sparse (1:2:2000, 9999, 1:1000, 1e4, 1e4);
  bm(i) = toc;
  pause (.010);
endfor


By putting everything in the next to last column this was clearly worst case, but it did beat the handcoded for loop.  Replacement code was


-          xcidx (0) = 0;
-          for (octave_idx_type j = 0; j < nc; j++)
-            xcidx (j+1) = j >= c(0) ? new_nz : 0;
+          std::fill_n (xcidx () + c(0) + 1, nc - c(0), new_nz);


I checked in the full cset here (http://hg.savannah.gnu.org/hgweb/octave/rev/584971932def).

Rik <rik5>
Group administrator
Mon 04 Sep 2017 08:38:08 AM UTC, comment #7: 

Rik, I'm not sure the gains will be great doing what you suggest. Creating an empty SparseRep and then resizing it once doesn't seem to be a large extra cost, as the creation of the empty SparseRep should be a low cost operation

D.

David Bateman <dbateman>
Group Member
Fri 01 Sep 2017 09:22:31 PM UTC, comment #6: 

Perhaps we should step away from this and look at it a different way.  Why does the constructor need to create an empty SparseRep and then use change_capacity to increase it?  I understand from a code clarity point-of-view that it guarantees a valid object is returned by the function regardless of which twists and turns the code takes through the various if/elseif branches.  However, given that this is a base library, I think it is more important to emphasize performance.

If we wait in the Sparse constructor until we know the number of elements, then we can use the this SparseRep constructor


    SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz = 0)
      : d (nz > 0 ? new T [nz] {} : nullptr),
        r (nz > 0 ? new octave_idx_type [nz] {} : nullptr),
        c (new octave_idx_type [nc+1] {}), nzmx (nz), nrows (nr),
        ncols (nc), count (1)
    { }


Because of the empty initializer list, the variables d, r, and c will all be filled with zeros.  That means that only the cidx entries from c(0) to nc need to be filled in rather than running over all columns nc.

Instead of the newly changed code from jwe:


  xcidx (0) = 0;
  for (octave_idx_type j = 0; j < nc; j++)
    xcidx (j+1) = j >= c(0) ? new_nz : 0;


It could be


// xcidx (0) = 0;  # no longer necessary, cidx is initialzed to all 0
std::fill_n (xcidx () + c(0) + 1, nc - c(0), new_nz);


This would also avoid some for loops and potential new[] operators in change_capacity.


Rik <rik5>
Group administrator
Wed 30 Aug 2017 06:56:36 PM UTC, comment #5: 

The currently applied patch has a for loop over 'nc' which is the number of columns which could be quite large for a sparse matrix, say 1e7.  David's suggestion would improve performance as the for loop would only be over the first c(0) blank columns which should be much less than 'nc' unless all of the data is on the right-hand side of the sparse matrix.

However, I couldn't get the code to work.  I tried this


      for (octave_idx_type i = 0; i <= c(0); i++)
        xcidx(i) = 0;


According to David, all values above should be 'nz' already from the call to change_capacity.  However, I used an assert and found that wasn't the case.


      assert (xcidx (c(0)+1) == std::max (nzm, new_nz));


change_capacity() is in Sparse.h


  void change_capacity (octave_idx_type nz)
  {
    if (nz < nnz ())
      make_unique (); // Unshare now because elements will be truncated.
    rep->change_length (nz);
  }


And change_length is in Sparse.cc


template <typename T>
void
Sparse<T>::SparseRep::change_length (octave_idx_type nz)
{
  for (octave_idx_type j = ncols; j > 0 && c[j] > nz; j--)
    c[j] = nz;


This would seem to set the last entries to the new_nz only on a truncation.  So, in the constructor which is growing the sparse array (nz = 0 to true nz0), maybe we have to keep the loop?  Comments?

Rik <rik5>
Group administrator
Tue 29 Aug 2017 01:59:47 PM UTC, comment #4: 

All values of xcidx will be initialised to nz by the call to change_capacity just before hand

D.

David Bateman <dbateman>
Group Member
Tue 29 Aug 2017 01:47:18 PM UTC, comment #3: 

I pushed this changeset:

http://hg.savannah.gnu.org/hgweb/octave/rev/520c37805969

This seems to work for me, but if there is a better fix, let me know.

John W. Eaton <jwe>
Group administrator
Tue 29 Aug 2017 01:34:12 PM UTC, comment #2: 

If I understand correctly, the xcidx vector has length nc+1 and with this code, if the specified column is less than nc, then the remaining elements would be uninitialized.  Is that OK?

Also, is there any check to ensure that c(0) > 0?  Maybe it can't be less than zero since it is supposed to be an element of an index vector?

Anyway, I came up with this instead:


  xcidx (0) = 0;
  for (octave_idx_type j = 0; j < nc; j++)
    xcidx (j+1) = j >= c(0) ? new_nz : 0;


but I don't know whether it is correct to set all the remaining xcidx elements to new_nz.

John W. Eaton <jwe>
Group administrator
Tue 29 Aug 2017 01:02:48 PM UTC, comment #1: 

I haven't done much with Octave for a long while, but as I wrote this code I had a quick look at it.

The specialisation for sparse column vectors was added to my code by Jaroslav, and he didn't take into account the possibility that the sparse column vector could in fact be a matrix. The fix is relatively easy. Looking at lines 362-363 of



http://hg.savannah.gnu.org/hgweb/octave/annotate/df49ac93f50c/liboctave/array/Sparse.cc


These implicitly assume that everything is in the first column. These should be replaced with


for (octave_idx_type i = 0; i < c (0) - 1; i++)
    xcidx (i) = 0;
xcidx (c(0) - 1) = new_nz;


I don't even have a octave source tree on my machine to test this so someone else will have to work this up as a hg chageset

D.

David Bateman <dbateman>
Group Member
Mon 28 Aug 2017 07:39:57 PM UTC, original submission:  

The following should produce a sparse matrix with non-zeros in the 2nd column.

sparse(1:2, 2, 1:2, 2, 2)


But it puts them in the first column.

> full(sparse(1:2, 2, 1:2, 2, 2))
ans =

   1   0
   2   0


In fact, the following will put the non-zeros in column 1 for any scalar value of j between 1 and 5.


sparse(1:5, j, 1:5, 5, 5)


So it appears that it is ignoring any scalar value passed in as a column index and just using 1 instead. A scalar row index, on the other hand, is handled correctly.

And Matlab also handles the scalar column index as expected.

Ray Zimmerman <rdzman>

 

(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 jwe (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by rdzman (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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-09-06 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2017-08-29 jwe StatusNone Ready For Test

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code