First non-zeros of a sparse array

Hi, I'm searching for the fastest way to get the index of the first non zero element of a sparse array. Only nonzeros have been "created" when I built this sparse matrix ==> searching the first nonzero is the same as searching the element with the smallest index.

The array are not very big (~5,000 maximum) but I have to process it millions of times. (find the first nonzero of each array)

Their length are equal, of the form n*(2^k) (n and k are constant), and I know that the first non-zero element is of the form m*(2^k) + 1. I thus tried to do it with a while loop, going through the array with (2^k) step, but it still takes time.

There is probably a way to find it faster.

Thanks!

Answer 1

[the cyclist]

This code takes about 2.5 seconds to run on my machine.

% Create a sparse array [10% filled]
x = sparse(5000,1);
x(randsample(5000,500)) = 1;

% Find the first element, a million times
idx = zeros(1000000,1);
tic
for n=1:1000000
    idx(n) = find(x,1);
end
toc

[Martin]

Thank you for your answers. 'find' seems to be effectively a great solution but it still does not work so well: the problem is that my data are not a list of arrays but a cell of 1000 sparse matrices:

P = cell(1,1000) % 1000 actions

And each P{i} is a sparse matrix, its size is about 5000*5000. (nState * nState)

I want to find, for each matrix, the first non-zero in each row, and store it in firstNonZero.

firstNonZero(nAction, nState) = 0;     % 1000 * 5000

for iAction = 1:nAction    % 1000 actions
  for iState = 1:nState    % 5000 states
    tempArray = P{iAction}(iState,:);   % Very slow !!!
    firstNonZero(iAction, iState) = find(tempArray,1); % Fast
       end
end

It's probably a stupid question again but it's now more about how to pick a row from a sparse matrix than taking its first nonzero. If I put the 2 relevant lines of the loop in one, it slows down.

[the cyclist]

You don't need to make the temp array:

firstNonZero(iAction, iState) = find(P{iAction}(iState,:),1);

Does that help?

[Martin]

No, unfortunately it does not help. Thanks anyway.

Answer 2

[per isakson]

Did you try

    ix = find( A ~= 0, 1, 'first' );

and something like

    ix = find( A("m*(2^k) + 1") ~= 0, 1, 'first' );

"non zero element" does that imply an integer array?

.

Doc says:

"[row,col] = find(X, ...) returns the row and column indices of the nonzero entries in the matrix X. This syntax is especially useful when working with sparse matrices."

[Richard Brown]

Or just find(A, 1)

Answer 3

[James Tursa]

You can use a mex routine for this, but you will need a C compiler to compile it. E.g., create a file called firstNonZeroMex.c with the contents below and then compile it with the mex function. Unfortunately, if you are on a 64-bit system then you will have to get and install a C compiler yourself since 64-bit MATLAB does not come with one installed. Not sure how this compares timing wise, but I would expect it to be fast compared to m-file methods that have to copy a lot of data to get at the rows. The mex function below does not copy any of the rows to get the job done.

/***********************************************************************
 * firstNonZeroMex returns first non-zero from each row of each sparse
 * matrix in a cell array of sparse matrices.
 *
 * Syntax:   F = firstNonZeroMex(C)
 *
 * Where:    C = Cell Array of sparse matrices, all the same row size
 *           F = First non-zero values of the C{i}
 *
 * If there is no non-zero in a particular row, returns 0 for that spot.
 * The result for C{i} row j will be in F(i,j).
 *
 * Programmer: James Tursa
 * Date: 25-July-2013
 ***********************************************************************/

#include "mex.h"
#ifndef MWSIZE_MAX
#  define mwSize int
#  define mwIndex int
#endif

void firstNonZero(double *pr, mxArray *cell, mwSize e, mwSize m);

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  mxArray *cell;
    mwSize i, m, n, e;
    double *pr;

  if( nrhs != 1 || !mxIsCell(prhs[0]) ) {
    mexErrMsgTxt("Need one cell array input");
  }
  if( nlhs > 1 ) {
    mexErrMsgTxt("Too many outputs");
  }
    e = mxGetNumberOfElements(prhs[0]);
  if( e == 0 ) {
    mexErrMsgTxt("Input is empty");
  }
  cell = mxGetCell(prhs[0],0);
  if( cell == NULL ) {
    mexErrMsgTxt("One or more of the cells is empty");
  }
  if( !mxIsSparse(cell) || !mxIsDouble(cell) ) {
    mexErrMsgTxt("One or more of the cells is not double sparse");
  }
  for( i=0; i<e; i++ ) {
    cell = mxGetCell(prhs[0],i);
      if( cell == NULL ) {
        mexErrMsgTxt("One or more of the cells is empty");
      }
      if( !mxIsSparse(cell) || !mxIsDouble(cell) ) {
        mexErrMsgTxt("One or more of the cells is not double sparse");
      }
    if( i == 0 ) {
          m = mxGetM(cell);
      plhs[0] = mxCreateDoubleMatrix(e,m,mxREAL);
      pr = mxGetPr(plhs[0]);
    } else {
      if( mxGetM(cell) != m ) {
        mexErrMsgTxt("One or more cells is mismatched in size from other cells");
      }
    }
    firstNonZero(pr, cell, e, m);
    pr++;
  }
}

void firstNonZero(double *pr, mxArray *cell, mwSize e, mwSize m)
{
    mwSize k, n, nrow;
    mwIndex *ir, *jc;
    mwIndex j, x, y;
    double *pc, *pj;

      n = mxGetN(cell);
      pc = mxGetData(cell);
      ir = mxGetIr(cell);
      jc = mxGetJc(cell);

  k = 0;
    for( y=0; y<n; y++ ) {
        nrow = jc[y+1] - jc[y];
        for( x=0; x<nrow; x++ ) {
      pj = pr + (*ir) * e;
      if( *pj == 0.0 ) {
        *pj = *pc;
        k++;
        if( k == m ) return;
      }
      ir++;
      pc++;
        }
    }
}

[Martin]

Thanks James. I'll try first to compute what I want at the moment I fill these matrices. If it is still slow, I'll try your solution. Thanks for your time.