First non-zeros of a sparse array
2 views (last 30 days)
Show older comments
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!
0 Comments
Accepted Answer
the cyclist
on 26 Jul 2013
Edited: the cyclist
on 26 Jul 2013
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
3 Comments
the cyclist
on 26 Jul 2013
You don't need to make the temp array:
firstNonZero(iAction, iState) = find(P{iAction}(iState,:),1);
Does that help?
More Answers (2)
per isakson
on 26 Jul 2013
Edited: per isakson
on 26 Jul 2013
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."
James Tursa
on 26 Jul 2013
Edited: James Tursa
on 26 Jul 2013
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++;
}
}
}
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!