How to quickly find the column index of the last non-zero element in all rows in a sparse matrix?

Hi, All,

I have a big sparse matrix A. I want to find out the column index of the Last non-zero element in all rows from the end. Here is my code:

reorderRow = [];

for jRow = 1 : length(A(:,1))

eee = find(A(jRow,:),1,'last');

reorderRow = [reorderRow eee];

end

For example, I have matrix A = [2 0 0;0 5 0;0 0 1;0 3 0;-1 0 -4]. My code gives the following result:

reorderRow = [1 2 3 2 3];

I am wondering if it is possible to obtain reorderRow without iterations. Thanks a lot.

Benson

Answer 1

[BobH]

reorderRow = arrayfun(@(X)  find(A(X,:),1,'last'), 1:size(A,1))
reorderRow =
     1     2     3     2     3

[John D'Errico]

As my other comments point out, IF the matrix is truly sparse, then the code you show can easily fail.

A = sprand(10000,10000,.0005);

reorderRow = arrayfun(@(X)  find(A(X,:),1,'last'), 1:size(A,1));
Error using arrayfun
Non-scalar in Uniform output, at index 166, output 1.
Set 'UniformOutput' to false. 

The failure arises since some rows can easily be empty of any elements. This can often be the case for a truly sparse matrix.

[Benson Gou]

Dear John,

Thanks a lot for your great help.

Benson

[Ameer Hamza]

Benson, I think you accepted the wrong answer. John's answer is under this answer.

[BobH]

Yes, Benson, please re-select John's answer

Answer 2

[John D'Errico]

It is important that if your matrix is sparse, that we use sparse matrix examples to test anything out.

Regardless, I would always do a comparison of other methods to a direct find. For example...

A = sprand(10000,10000,.0005);

So, on average, 5 elements per row. If some rows have no non-zeros, then this will return zero as the index of the last element.

tic
[ir,ic] = find(A);
maxcolind = accumarray(ir,ic,[10000,1],@max);
toc
Elapsed time is 0.026327 seconds.

This will not fail, even if some rows are empty of non-zeros.

sum(maxcolind == 0)
ans =
    59

The random example I created had 59 rows that were entirely zero.

If you want to find the array reorderRow as you asked for, then just do this:

r = find(maxcolind);
reorderRow = [r,maxcolind(r)];

[Benson Gou]

Dear John,

Thanks a lot. Your code took 0.002112 s for a sparse matrix 974 x 632. It works well.

Benson

Answer 3

[James Tursa]

You can use a mex routine for this and avoid all temporary data copies. E.g., a direct approach:

/*   File:  find_last_nz.c
 *
 *   C = find_last_nz(S)
 *
 *   Finds column number of last non-zero element in each row of sparse matrix
 *   S = sparse matrix (double or logical)
 *   C = full column vector containing the columns numbers for each row
 *
 *   Programmer: James Tursa
 */

/* Includes ----------------------------------------------------------- */

#include "mex.h"

/* Gateway ------------------------------------------------------------ */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *Cpr;
    mwIndex *Ir, *Jc;
    size_t i, j, imax, m, n;
    
    if( nlhs > 1 ) {
        mexErrMsgTxt("Too many outputs.");
    }
    if( nrhs != 1 || !mxIsSparse(prhs[0]) ) {
        mexErrMsgTxt("Need exactly one sparse input.");
    }
    Ir = mxGetIr(prhs[0]);
    Jc = mxGetJc(prhs[0]);
    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);
    plhs[0] = mxCreateDoubleMatrix(m,1,mxREAL);
    Cpr = (double *) mxGetData(plhs[0]);
    for( j=0; j<n; j++ ) {
        imax = Jc[j+1] - Jc[j];
        for( i=0; i<imax; i++ ) {
            Cpr[*Ir++] = j+1;
        }
    }
}

And a sample run:

>> mex find_last_nz.c
Building with 'Microsoft Windows SDK 7.1 (C)'.
MEX completed successfully.
>> A = rand(10).*(rand(10)<.2)
A =
    0.6443         0    0.3111         0         0         0    0.0377         0         0         0
         0         0         0         0         0         0         0         0         0         0
    0.8116         0         0    0.6028    0.8010         0         0    0.4942         0         0
    0.5328         0         0         0         0         0         0         0         0         0
         0         0    0.9049    0.2217         0         0    0.0987         0         0         0
    0.9390         0         0    0.1174         0         0         0    0.9037         0    0.1679
         0         0         0         0         0    0.6791         0         0         0         0
         0         0         0    0.3188         0    0.3955         0         0         0         0
         0    0.2277         0         0         0         0    0.1366         0         0    0.5005
    0.5870    0.4357         0         0         0    0.9880         0         0    0.5767    0.4711
>> find_last_nz(sparse(A))
ans =
     7
     0
     8
     1
     7
    10
     6
     6
    10
    10

What is presented above is a one-pass approach that only looks at the internal sparse indexing arrays.

It is possible that it might be faster to transpose the matrix first so that I could pick off the column numbers directly from the internal arrays, but that would require a data copy and I haven't written and tested this approach. That is, picking off the last non-zero row number for each column would be an extremely fast operation in a mex routine, but for your case the matrix would have to be transposed first which would probably dominate the run-time so I haven't bothered to write and test it.

EDIT

A timing comparison with the m-code method posted by John:

>> A = sprand(10000,10000,.001);
>> 
>> tic;[ir,ic] = find(A);maxcolind = accumarray(ir,ic,[10000,1],@max);toc
Elapsed time is 0.029285 seconds.
>> tic;[ir,ic] = find(A);maxcolind = accumarray(ir,ic,[10000,1],@max);toc
Elapsed time is 0.030831 seconds.
>> tic;[ir,ic] = find(A);maxcolind = accumarray(ir,ic,[10000,1],@max);toc
Elapsed time is 0.028965 seconds.
>> 
>> tic;C=find_last_nz(A);toc
Elapsed time is 0.000260 seconds.
>> tic;C=find_last_nz(A);toc
Elapsed time is 0.000253 seconds.
>> tic;C=find_last_nz(A);toc
Elapsed time is 0.000278 seconds.
>> 
>> isequal(maxcolind,C)
ans =
  logical
   1

So, the mex routine is about 100x faster for this particular example.

[Benson Gou]

Hi, James,

Thanks for your great help. But I do not know how to connect your c code with my Matlab code.

Benson

[Matt J]

Just do,

>> mex find_last_nz.c

[Benson Gou]

Dear James,

Thanks a lot for your great. It seems that your code can reach a very small running time. Now I am learning how to call it in my Matlab code.

Thanks a lot again.

Benson

[Ameer Hamza]

Benson, you just need to copy the James's code and paste it in a file named find_last_nz.c. Then run

mex find_last_nz.c

After that, you can call this function, like any other MATLAB function

find_last_nz(A)

[Matt J]

Benson, you just need to copy the Matt's code and paste it in a file named find_last_nz.c. Then run

James's code, I think you mean.

[Benson Gou]

Dear Matt,

Good morning!

Thanks a lot for your great help. You have a good weekend!

Benson

Answer 4

[Ameer Hamza]

For starters, do preallocation, It can speed up the execution if the matrix A has many rows. The following code compares the execution time for your code, BobH's code, and the code with pre-allocation. I intentionally made the A matrix with 10000 rows to show the difference

A = randi([1 5], 10000, 3);

% your code
tic
reorderRow = [];
for jRow = 1 : length(A(:,1))
    eee = find(A(jRow,:),1,'last');
    reorderRow = [reorderRow eee];
end
toc

% arrayfun method
tic
reorderRow = arrayfun(@(X)  find(A(X,:),1,'last'), 1:length(A));
toc

% pre-allocating the array
tic
reorderRow = zeros(1, size(A,1));
for jRow = 1 : length(A(:,1))
    eee = find(A(jRow,:),1,'last');
    reorderRow(jRow) = eee;
end
toc

Result

Elapsed time is 0.068404 seconds. % your code (slowest)
Elapsed time is 0.048778 seconds. % array fun (medium)
Elapsed time is 0.015521 seconds. % pre-allocation (fastest)

Pre-allocation is about 3.5 times faster than the original code.

[John D'Errico]

Should I point out that your example matrix is not sparse? As such, it could lead to very different timings, when used on truly sparse matrices. The point being that any operations on sparse matrices tand to have different timings than they do on fiull matrices. Even something as simple as extracting the i'th row of a sparse matrix can be a relatively costly thing, since the elements in a sparse matrix are stored in the wrong sequence, going down the columns.

[Ameer Hamza]

John, Thanks for pointing out. I noticed that before, but all of these routines are essentially calling the same functions find(A(jRow,:),1,'last') with similar options, albeit differently. So, I thought this should be a fair comparison. If there is some other method that does not involve find(), than sparsity might play a role in comparison.

I tried with the sparse matrix as given by OP, and results appear to be similar.

A = [2 0 0;0 5 0;0 0 1;0 3 0;-1 0 -4];
A = repmat(A, 2000, 1); % <--- simply repeated the matrix given by OP so have same sparsity

% your code
tic
reorderRow = [];
for jRow = 1 : length(A(:,1))
    eee = find(A(jRow,:),1,'last');
    reorderRow = [reorderRow eee];
end
toc

% arrayfun method
tic
reorderRow = arrayfun(@(X)  find(A(X,:),1,'last'), 1:length(A));
toc

% pre-allocating the array
tic
reorderRow = zeros(1, size(A,1));
for jRow = 1 : length(A(:,1))
    eee = find(A(jRow,:),1,'last');
    reorderRow(jRow) = eee;
end
toc

Result:

Elapsed time is 0.076666 seconds.
Elapsed time is 0.056438 seconds.
Elapsed time is 0.015933 seconds.

[John D'Errico]

You also need to watch out for empty rows, as that will cause your code to fail, with this error:

"Unable to perform assignment because the left and right sides have a different number of elements."

And if the matrix is indeed sparse, that is an entirely possible event.

[Ameer Hamza]

Yes, I saw your solution. That's indeed brilliant.

[Benson Gou]

Dear Ameer,

Thanks a lot for your great help.

Benson

Answer 5

[Matt J]

[~,idx]=max(fliplr(logical(A)),[],2);
result=size(A,2)+1-idx;

[Matt J]

And possibly you might add,

result(idx==0)=nan; %The row has no non-zeros

[Benson Gou]

Dear Matt,

Thanks a lot. I compared your code for the "first non-zero element for columns" and the "last non-zero element for rows", I found the later takesa little more time (0.000583 s vs 0.002868 s) for 974 x 632 sparse matrix. But it works well. Thanks.

Benson

[Matt J]

If the difference really matters, you might consider transposing the matrix and just working on columns exclusively.

[Benson Gou]

Hi, Matt,

I tried to use the transpose of matrix A and apply your code. But after the transposing, we need to find the index of the last non-zero element in each column instead of the first non-zero element. But logical only returns the index of the first non-zero element of each column.

The time does matter because we have to run it for more than 1000 times.

Thanks.

Benson

[Matt J]

You need to use flipud() so that "first" becomes "last", similar to how fliplr() did this in my solution above for row-wise search.

[Benson Gou]

Thanks, Matt. It is very helpful. I will try your idea.

Best regards,

Benson