# bug #53198: nnz could be improved for diagonal matrices

 Submitter: Rik Submitted: Wed 21 Feb 2018 06:11:33 AM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Performance Status: Fixed Assigned to: None Originator Name: Open/Closed:  * Closed Release:  * dev Operating System:  * Any Fixed Release: None Planned Release: None
* Mandatory Fields

Add a New Comment Rich Markup

## Discussion

Mon 18 Jun 2018 07:53:11 PM UTC, comment #9:

I created a few bug reports to address the issues I've noted, bug #54142, bug #54143, bug #54144.  There could be another one or two as well, but let's discontinue this thread and continue in the individual new reports since the other issues aren't related to nnz.

Dan Sebald <sebald>
Mon 18 Jun 2018 06:45:50 PM UTC, comment #8:

Oh, and I forgot to mention one of the first things that caught my eye, and why I started the examples with a sparse matrix.  In the issorted() documentation it states

This function does not support sparse matrices.

Again, sorting a sparse matrix is simply a case of changing the sparse-matrix's indeces appropriately.  In fact, it's already implemented:

octave:33> sort(X)
ans =

Compressed Column Sparse (rows = 3, cols = 4, nnz = 2 [17%])

(3, 1) ->  4
(3, 2) ->  5

So why can't issorted() work for sparse matrices?

Dan Sebald <sebald>
Mon 18 Jun 2018 06:29:48 PM UTC, comment #7:

I'll give it a try.  I'm a bit busy this summer, so it may take a while even though it's not that big of a changeset.

I'm a little confused with the current issorted() for full matrices.  The following all looks fine:

octave:20> I = [1 1 2]; J = [1 1 2]; SV = [3 4 5];
octave:21> X = sparse (I, J, SV, 3, 4, "unique")
X =

Compressed Column Sparse (rows = 3, cols = 4, nnz = 2 [17%])

(1, 1) ->  4
(2, 2) ->  5

octave:22> full(X)
ans =

4   0   0   0
0   5   0   0
0   0   0   0

octave:23> sort(full(X))
ans =

0   0   0   0
0   0   0   0
4   5   0   0

octave:24> sort(full(X),2)
ans =

0   0   0   4
0   0   0   5
0   0   0   0

But why can't issorted() act on a matrix?

octave:25> issorted(sort(full(X)))

I would think that issorted() is a counterpart to sort(), so it should be able to act on something output by sort().

Furthermore, there is this "rows" option, which refers to sortrows().  Now, I understand the idea of a sortrows().  It's a vector-based sort that maintains order across rows; a bit like a spreadsheet sort, if you follow.  But it seems to me that in order to test whether this is a valid output of sortrows(), one needs to know how the sorting of sortrows() was done, i.e., the "C" of sortrows() documentation.  Also, note that this works

octave:30> issorted(sort(full(X)),"rows")
ans = 1

yet without the "rows" it's the error above.  And what about a "cols" counterpart?  Why should vector-sorting be allowed only in the "rows" dimension?

Dan Sebald <sebald>
Mon 18 Jun 2018 03:50:39 PM UTC, comment #6:

Using the 80/20 rule, I think the biggest gains would be made by changing sort and issorted.  I'm okay with leaving reshape and permute as they are.

There is some question whether this should be done in the interpreter library or in liboctave but we could punt on that for now and just do it in the interpreter.

Do you want to go ahead and convert the pseudocode to a patch?

Rik <rik5>
Sun 17 Jun 2018 06:46:45 PM UTC, comment #5:

That's nice.

I'm looking back at the comments I made and remembered some of this.  There is the issue that other uses of to_dense() might be avoided.  That is, rather than use to_dense() and then the slower non-sparse methods, construct a dense() zeros matrix and fill in only the necessary fields.  Most of these won't be major speed improvements like nnz, and I wonder how often these others get used.  Nonetheless, here are some comments about each (actually making the mods means altering a header file which takes considerable time to recompile):

RESHAPE

octave_value reshape (const dim_vector& new_dims) const

This wouldn't be a huge speed up because I would think reshape could be somewhat fast, O(N) (I haven't tried it).  However, what this routine is doing ultimately is just distributing those diagonal values to a particular set of locations that could probably be handled with a simple 2D algorithm that computes the indices appropriately.  That is,

1) Determine the ultimate shape of the output matrix.
2) Construct an all-zeros matrix of that size.
3) Take the diagonal vector and distribute it within that matrix, something like

x_out = zeros(N_i,N_j);
for (i=0; i<N_i; i+=S_i)
{
for (j=0; j<N_j; j+=S_j)
{
x_out(i,j) = x_diag(i_diag++);
}
}

PERMUTE

octave_value permute (const Array<int>& vec, bool inv = false) const
{
if (vec.numel () == 2
&& ((vec.xelem (0) == 1 && vec.xelem (1) == 0)
|| (vec.xelem (0) == 0 && vec.xelem (1) == 1)))
return DMT (matrix);
else
}

This one probably doesn't matter too much.  Even in the cases where the number of dimensions are expanded, e.g., permute (X, [3, 1, 2]), where X is a 2D diagonal matrix, the number of elements remain the same.  So, an indexing formula approach could be used here as well, but to_dense() may be fine.

SORT

octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
sortmode mode = ASCENDING) const

This is the one that originally caught my attention.  The algorithm is simple.

1) Construct an all-zeros matrix of the same size.
2) Take the diagonal vector and distribute it either to the starting or ending position of the colum (or row) of the zeros matrix.  E.g.,

x_out = zeros(N_i,N_j);
for (ivec=0; ivec<N_vec; ivec++)
{
if (xsparse(ivec) < 0)
{
if (col)
x_out(0, ivec) = x_diag(ivec);
else
x_out(ivec, 0) = x_diag(ivec);
}
else
{
if (col)
x_out(N_i-1, ivec) = x_diag(ivec);
else
x_out(ivec, N_j-1) = x_diag(ivec);
}
}

ISSORTED

sortmode issorted (sortmode mode = UNSORTED) const

This one could be a considerable savings.  Isn't the answer to this simply dependent on the nature of the diagonal values given we already know that everything else is zero?  There are a few exceptions:

1) If the matrix size is 1 x N, the matrix is sorted along columns.
2) If the matrix size is M x 1, the matrix is sorted along rows.
3) If the matrix size is 2 x 2 and xdiag(0) < 0 and xdiag(1) > 0, the matrix is sorted.
3) If the diagonal vector is all zeros, the matrix is sorted.
4) Otherwise, the matrix is not sorted.

Dan Sebald <sebald>
Sun 17 Jun 2018 04:19:26 PM UTC, comment #4:

I pushed the modifide patch here (http://hg.savannah.gnu.org/hgweb/octave/rev/920299ced721).  Dan was correct that there is no need to check for whether the object is a diagonal matrix inside the diagonal class.  Speed up is 1500X for the test script.

Rik <rik5>
Wed 21 Feb 2018 09:25:54 PM UTC, comment #3:

I've applied a second test with

N = 1e2;

f = zeros(N,1);
for i=1:N
s = diag(1,1,5000);
tic;
z = nnz(s);
f(i)=toc;
endfor

z
sum(f)

And I don't see any consequence of the patch

octave:1> % Pre-patch
octave:1> test
ans =  17.969
octave:2> test2
z =  1
ans =  0.0031841

octave:1> % Post-patch
octave:1> test
ans =  0.0032768
octave:2> test2
z =  1
ans =  0.0043147

I take it that in this scenario:

if (k == 0)
// Returns Diag2Array<T> with nnz <= 1.
retval = matrix.build_diag_matrix ();
else

the matrix representation is some other class that has a fast implementation of nnz().  So that's good.

One other note: be sure to check changesets for any dangling whitespace.  If you are running linux, type "hg diff" at a bash shell and you are likely to see highlighted red blocks wherever there is whitespace after the last character.

Dan Sebald <sebald>
Wed 21 Feb 2018 08:01:14 PM UTC, comment #2:

I've looked at the patch.  The conditional test

if (is_diag_matrix ())

in the patch appears to be superfluous.  Isn't the class, by definition, diagonal?  Wouldn't simply extracting the diagonal, i.e., collapsing to a vector, be sufficient?  I.e.,

octave_idx_type nnz (void) const
{
return diag ().nnz ();
}

However, there is some special code in the ::diag() method:

if (matrix.rows () == 1 || matrix.cols () == 1)
{
// Rather odd special case.  This is a row or column vector
// represented as a diagonal matrix with a single nonzero entry, but
// Fdiag semantics are to product a diagonal matrix for vector
// inputs.
if (k == 0)
// Returns Diag2Array<T> with nnz <= 1.
retval = matrix.build_diag_matrix ();
else
// Returns Array<T> matrix
retval = matrix.array_value ().diag (k);
}
else
// Returns Array<T> vector
retval = matrix.extract_diag (k);
return retval;

Is that special case something to be concerned about?  I think that would arise for the following sort of construct:

octave:17> a = diag(1,1,10000);
octave:18> isdiag(a)
ans = 1
octave:19> size(a)
ans =

1   10000

octave:20> size(diag(a,0))
ans =

10000   10000

(Note how the application of diag() with k=0, the default in calling the ::diag() member function, expands the matrix in this case, not contracts it.)

Also, throughout ov-base-diag.h are several to_dense() calls.  Maybe they all could be addressed with one patch here.  For example, there is the method "sort" that uses to_dense() first, but that too could be improved speed-wise noting that the operation boils down to pushing the diagonal value either to the first or last row, WLOG, e.g.,

octave:10> sort(diag([1 -2 3]))
ans =

0  -2   0
0   0   0
1   0   3

(if it is sort-by-row, then change the rule for the other direction).  One would just start with the appropriately sized zeros matrix and then index through the diagonal vector sending the values to the top or bottom row based on one less-than-zero comparison.

Dan Sebald <sebald>
Wed 21 Feb 2018 06:03:05 PM UTC, comment #1:

Indeed the time difference is huge.

I tried the following test script:

N = 1e2;

f = zeros(N,1);
for i=1:N
r = rand(5000,1);
s = diag(r);
tic;
z = nnz(s);
f(i)=toc;
endfor

sum(f)

Got a total time of 0.0039184 seconds with the patch applied and 11.160 seconds without the patch.

Hope this helps :)

(file #43379, file #43380)

Sahil <batterylow>
Wed 21 Feb 2018 06:11:33 AM UTC, original submission:

nnz for diagonal matrices is implemented in libinterp/octave-value/ov-base-diag.h.  The code is simple

However, this isn't very efficient.  It converts and then caches the full version of the matrix and then calls nnz on that full matrix.  Besides the conversion, nnz on a full matrix means looking at every element.  For a diagonal matrix nnz should just need to look at the elements on the diagonal.

For a square matrix the difference is O(N^2) for full matrix versus O(N) for just the diagonal.

Rik <rik5>

## Attached Files

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

Attach Files:

Comment:

Attached Files
file #43379:  53198.diff added by batterylow (909B - text/x-patch)
file #43380:  test.m added by batterylow (118B - text/x-objcsrc)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by sebald (Posted a comment)
• -email is unavailable- added by batterylow (Updated the item)
• -email is unavailable- added by rik5 (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.