Thread Subject: speeding up: min(A(A>0))

Hello all,

I have a while minA<threshold loop and inside I minA is
defined as follows:

minA=min(A(A>0));

where A is a positive matrix (aprox size 80 by 40) that
has several zeros. The number of zeros is updated in the
loop.
Using the profiler I realized that I am computing minA
about 16.5e6 times and it takes about 1300 secs in my pc
which is about 50% of the total running time of my script.
Is there any way to speed up this code?

Thanks!!

"carmelo " <km.nospam.scriba.nospam@yahoo.esnospam> writes:

> Hello all,
>
> I have a while minA<threshold loop and inside I minA is
> defined as follows:
>
> minA=min(A(A>0));
>
> where A is a positive matrix (aprox size 80 by 40) that
> has several zeros. The number of zeros is updated in the
> loop.
> Using the profiler I realized that I am computing minA
> about 16.5e6 times and it takes about 1300 secs in my pc
> which is about 50% of the total running time of my script.
> Is there any way to speed up this code?

In isolation, no. There might be a way to avoid searching the full
array for zeros by changing the algorithm, for instance by computing
only the A elements that change.

If you always set A elements to zero, that would be equivalent to
removing those elements from a list of nonzero elements that would get
smaller and smaller each time.

You could short-circuit the minimum by comparing only the changed values
against the current minimum. If changed A values are all greater than
the current minimum, no need to recompute.

If you always add A elements, then simply keep a running minimum.

The possibilities are endless... the specifics depend on your algorithm.

-Peter

> > Hello all,
> >
> > I have a while minA<threshold loop and inside I minA
is
> > defined as follows:
> >
> > minA=min(A(A>0));
> >
> > where A is a positive matrix (aprox size 80 by 40)
that
> > has several zeros. The number of zeros is updated in
the
> > loop.
> > Using the profiler I realized that I am computing minA
> > about 16.5e6 times and it takes about 1300 secs in my
pc
> > which is about 50% of the total running time of my
script.
> > Is there any way to speed up this code?
>
> In isolation, no. There might be a way to avoid
searching the full
> array for zeros by changing the algorithm, for instance
by computing
> only the A elements that change.
>
> If you always set A elements to zero, that would be
equivalent to
> removing those elements from a list of nonzero elements
that would get
> smaller and smaller each time.
>
> You could short-circuit the minimum by comparing only
the changed values
> against the current minimum. If changed A values are
all greater than
> the current minimum, no need to recompute.
>
> If you always add A elements, then simply keep a running
minimum.
>
> The possibilities are endless... the specifics depend on
your algorithm.
>
> -Peter


Thanks Peter,

here is the algorithm

minA=min(A(A>0));
while minprop<threshold
    [row,col]=find(A==minprop,1);
    nameclasses_row=nameclasses(A(row,:)~=0);
    for iiclass=1:numberclasses-1
        candidate_replacer=ordering(nameclasses
(col),iiclass+1);
        if any(candidate_replacer==nameclasses_row)
            % find col of candidate_replacer
            col_replacer=find
(nameclasses==candidate_replacer);
            % do replacement
            A(row,col_replacer)=A(row,col_replacer)+A
(row,col);
            A(row,col)=0;
            break
        end
    end
    minA=min(A(A>0));
end

I search for the min element in A, then i use the
matrix "ordering" to search for the column (in the row
where minA was found) where minA should be added. the
matrix "ordering" gives a secuential order of cols to add
minA. This is repeated until all the elements in A are
above a certain threshold.

so i guess that the solution is to only update the indices
of A>0 if A has changed....

"carmelo " <km.nospam.scriba.nospam@yahoo.esnospam> wrote
in message <fqm7k2$nug$1@fred.mathworks.com>...
> Hello all,
>
> I have a while minA<threshold loop and inside I minA is
> defined as follows:
>
> minA=min(A(A>0));
>
> where A is a positive matrix (aprox size 80 by 40) that
> has several zeros. The number of zeros is updated in the
> loop.
> Using the profiler I realized that I am computing minA
> about 16.5e6 times and it takes about 1300 secs in my pc
> which is about 50% of the total running time of my
script.
> Is there any way to speed up this code?
>
> Thanks!!

If finding non-zero indices is your only concern, you may
find (~~A) is quicker than (A>0). On my old PC running
2006b the (~~A) approach was about twice the speed of (A>0).

carmelo <km.nospam.scriba.nospam@yahoo.esnospam> wrote:
> Hello all,
>
> I have a while minA<threshold loop and inside I minA is
> defined as follows:
>
> minA=min(A(A>0));
>
> where A is a positive matrix (aprox size 80 by 40) that
> has several zeros.

What you are doing is already pretty good, and I don't see much scope for
improvement without resorting to some C code. That said, the following
does seem to help on examples I have tried:

BEGIN-----CUT HERE-----
/*
 * minnz.c - Find minimum non-zero value in an array
 *
 * function m = minnz(x)
 *
 * Tristram Scott
 * 05/03/2008
 *
 * $Id: minnz.c,v 1.2 2008/03/05 16:57:38 tristram Exp $
 */
#include <stdlib.h>
#include <math.h>
#include "mex.h"
/* the gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] ) {
const mxArray *tmp1;
double m, *mout, *x;
int i, n;

if(nrhs!=1) {
mexErrMsgTxt("One input required.");
}

tmp1 = prhs[0];
n = mxGetM(tmp1) * mxGetN(tmp1);
x = mxGetPr(tmp1);
m = 1.0/0;
for (i=0;i<n;i++) {
if ((x[i] < m) && (x[i])) {
m = x[i];
}
}

plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
mout = mxGetPr(plhs[0]);
*mout = m;

}
END-----CUT HERE-----




--
Dr Tristram J. Scott
Energy Consultant

tristram.scott@ntlworld.com (Tristram Scott) wrote in
message <zjAzj.63245$os2.35954@newsfe3-win.ntli.net>...
> carmelo <km.nospam.scriba.nospam@yahoo.esnospam> wrote:
> > Hello all,
> >
> > I have a while minA<threshold loop and inside I minA is
> > defined as follows:
> >
> > minA=min(A(A>0));
> >
> > where A is a positive matrix (aprox size 80 by 40) that
> > has several zeros.
>
> What you are doing is already pretty good, and I don't
see much scope for
> improvement without resorting to some C code. That said,
the following
> does seem to help on examples I have tried:
>
> BEGIN-----CUT HERE-----
> /*
> * minnz.c - Find minimum non-zero value in an array
> *
> * function m = minnz(x)
> *
> * Tristram Scott
> * 05/03/2008
> *
> * $Id: minnz.c,v 1.2 2008/03/05 16:57:38 tristram Exp $
> */
> #include <stdlib.h>
> #include <math.h>
> #include "mex.h"
> /* the gateway function */
> void mexFunction( int nlhs, mxArray *plhs[],
> int nrhs, const mxArray *prhs[] ) {
> const mxArray *tmp1;
> double m, *mout, *x;
> int i, n;
>
> if(nrhs!=1) {
> mexErrMsgTxt("One input required.");
> }
>
> tmp1 = prhs[0];
> n = mxGetM(tmp1) * mxGetN(tmp1);
> x = mxGetPr(tmp1);
> m = 1.0/0;
> for (i=0;i<n;i++) {
> if ((x[i] < m) && (x[i])) {
> m = x[i];
> }
> }
>
> plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
> mout = mxGetPr(plhs[0]);
> *mout = m;
>
> }
> END-----CUT HERE-----
>
>
>
>
> --
> Dr Tristram J. Scott
> Energy Consultant

Alternativley, and a bit more tersely:


i = mxGetM(tmp1)*mxGetN(tmp1) -1;
m = x[i];
while(i--) if ((x[i] < m) && (x[i])) m = x[i];

"Steve Amphlett" <Firstname.Lastname@Where-I-Work.com>
wrote in message <fqml72$reg$1@fred.mathworks.com>...
> tristram.scott@ntlworld.com (Tristram Scott) wrote in
> message <zjAzj.63245$os2.35954@newsfe3-win.ntli.net>...
> > carmelo <km.nospam.scriba.nospam@yahoo.esnospam>
wrote:
> > > Hello all,
> > >
> > > I have a while minA<threshold loop and inside I minA
is
> > > defined as follows:
> > >
> > > minA=min(A(A>0));
> > >
> > > where A is a positive matrix (aprox size 80 by 40)
that
> > > has several zeros.
> >
> > What you are doing is already pretty good, and I don't
> see much scope for
> > improvement without resorting to some C code. That
said,
> the following
> > does seem to help on examples I have tried:
> >
> > BEGIN-----CUT HERE-----
> > /*
> > * minnz.c - Find minimum non-zero value in an array
> > *
> > * function m = minnz(x)
> > *
> > * Tristram Scott
> > * 05/03/2008
> > *
> > * $Id: minnz.c,v 1.2 2008/03/05 16:57:38 tristram Exp
$
> > */
> > #include <stdlib.h>
> > #include <math.h>
> > #include "mex.h"
> > /* the gateway function */
> > void mexFunction( int nlhs, mxArray *plhs[],
> > int nrhs, const mxArray *prhs[] ) {
> > const mxArray *tmp1;
> > double m, *mout, *x;
> > int i, n;
> >
> > if(nrhs!=1) {
> > mexErrMsgTxt("One input required.");
> > }
> >
> > tmp1 = prhs[0];
> > n = mxGetM(tmp1) * mxGetN(tmp1);
> > x = mxGetPr(tmp1);
> > m = 1.0/0;
> > for (i=0;i<n;i++) {
> > if ((x[i] < m) && (x[i])) {
> > m = x[i];
> > }
> > }
> >
> > plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
> > mout = mxGetPr(plhs[0]);
> > *mout = m;
> >
> > }
> > END-----CUT HERE-----
> >
> >
> >
> >
> > --
> > Dr Tristram J. Scott
> > Energy Consultant
>
> Alternativley, and a bit more tersely:
>
>
> i = mxGetM(tmp1)*mxGetN(tmp1) -1;
> m = x[i];
> while(i--) if ((x[i] < m) && (x[i])) m = x[i];
>

Thanks Tristam, Thanks Peter,

I do not know much about C and i have never created MEX
files.
i put the code between the BEGIN and END keywords in a
text file. i saved it as minnz.c and i run mex minnz.c
then i got the following error msg:

lcc preprocessor warning: minnz.c:39 No newline at end of
file
Error minnz.c: 28 Division by zero
1 errors, 1 warnings

could you please help me with this?

"Tristram Scott" <tristram.scott@ntlworld.com> wrote in message
news:zjAzj.63245$os2.35954@newsfe3-win.ntli.net...
> carmelo <km.nospam.scriba.nospam@yahoo.esnospam> wrote:
>> Hello all,
>>
>> I have a while minA<threshold loop and inside I minA is
>> defined as follows:
>>
>> minA=min(A(A>0));
>>
>> where A is a positive matrix (aprox size 80 by 40) that
>> has several zeros.
>
> What you are doing is already pretty good, and I don't see much scope for
> improvement without resorting to some C code. That said, the following
> does seem to help on examples I have tried:

*snip MEX code*

That works for most cases. Being somewhat nitpicky, though, I see that
minnz(NaN) returns Inf and minnz([]) and minnz(0) also return Inf. NaN is
not equal to 0, so the code should return NaN in the first case, and in the
latter cases the result should be [], as the matrices for those cases have
no nonzero elements.

--
Steve Lord
slord@mathworks.com

Steven Lord <slord@mathworks.com> wrote:
>
> "Tristram Scott" <tristram.scott@ntlworld.com> wrote in message
> news:zjAzj.63245$os2.35954@newsfe3-win.ntli.net...
>> carmelo <km.nospam.scriba.nospam@yahoo.esnospam> wrote:
>>> Hello all,
>>>
>>> I have a while minA<threshold loop and inside I minA is
>>> defined as follows:
>>>
>>> minA=min(A(A>0));
>>>
>>> where A is a positive matrix (aprox size 80 by 40) that
>>> has several zeros.
>>
>> What you are doing is already pretty good, and I don't see much scope for
>> improvement without resorting to some C code. That said, the following
>> does seem to help on examples I have tried:
>
> *snip MEX code*
>
> That works for most cases. Being somewhat nitpicky, though, I see that
> minnz(NaN) returns Inf and minnz([]) and minnz(0) also return Inf. NaN is
> not equal to 0, so the code should return NaN in the first case, and in the
> latter cases the result should be [], as the matrices for those cases have
> no nonzero elements.
>

Well Steve, it's nice to know that someone reads my rubbish from time to
time. Certainly the code should check for those trivial cases. Perhaps
I'll tackle that in the morning...

Thanks for your comments. It is worth the effort to make routines like
this more robust.

--
Dr Tristram J. Scott
Energy Consultant

carmelo <km.nospam.scriba.nospam@yahoo.esnospam> wrote:
> I do not know much about C and i have never created MEX
> files.
> i put the code between the BEGIN and END keywords in a
> text file. i saved it as minnz.c and i run mex minnz.c
> then i got the following error msg:
>
> lcc preprocessor warning: minnz.c:39 No newline at end of
> file

Put a blank line at the end.

> Error minnz.c: 28 Division by zero
> 1 errors, 1 warnings

I can never remember what infinity is called, so I had it calculate one for
me.

A bit of digging about gives me the following:

m = INFINITY;

This seems to be correct under Solaris, but I don't know what lcc is going
to say about it.

--
Dr Tristram J. Scott
Energy Consultant

Steve Amphlett <Firstname.Lastname@where-i-work.com> wrote:
>
> Alternativley, and a bit more tersely:
>
>
> i = mxGetM(tmp1)*mxGetN(tmp1) -1;
> m = x[i];
> while(i--) if ((x[i] < m) && (x[i])) m = x[i];
>

That might just go a bit quicker if the test of (i--) is non-zero is
quicker than the test of (i<n). Good idea. Also, though, I'd be inclined
to use m = INFINITY, in case x[end] is zero. And see comments elsewhere in
this thread on other cases to test for.

--
Dr Tristram J. Scott
Energy Consultant

Tristram Scott <tristram.scott@ntlworld.com> wrote:
> Steven Lord <slord@mathworks.com> wrote:
>>
>> "Tristram Scott" <tristram.scott@ntlworld.com> wrote in message
>> news:zjAzj.63245$os2.35954@newsfe3-win.ntli.net...

[snip]

>>
>> That works for most cases. Being somewhat nitpicky, though, I see that
>> minnz(NaN) returns Inf and minnz([]) and minnz(0) also return Inf. NaN is
>> not equal to 0, so the code should return NaN in the first case, and in
>> the
>> latter cases the result should be [], as the matrices for those cases have
>> no nonzero elements.
>>
>
> Well Steve, it's nice to know that someone reads my rubbish from time to
> time. Certainly the code should check for those trivial cases. Perhaps
> I'll tackle that in the morning...
>

Well, the code has grown, but with little loss of efficiency in the tests I
have done.

BEGIN---cut here---
/*
 * minnz.c - Find minimum non-zero value in an array
 *
 * function m = minnz(x)
 *
 * Tristram Scott
 * 05/03/2008
 *
 * $Id: minnz.c,v 1.3 2008/03/06 09:02:55 tristram Exp $
 */
#include <stdlib.h>
#include <math.h>
#include "mex.h"
/* the gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] ) {
const mxArray *tmp1;
double m, *mout, *x;
int f, n;

if(nrhs!=1) {
mexErrMsgTxt("One input required.");
}

tmp1 = prhs[0];
n = mxGetM(tmp1) * mxGetN(tmp1);
if (!n) {
plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
return;
}
x = mxGetPr(tmp1);
f = 0; /* Flag */
while (n--) {
if (x[n]) {
m = x[n];
f = 1; /* Found a non-zero value */
if (!isnan(x[n])) {
break;
}
}
}
if (n>=0) {
while (n--) {
if ((x[n] < m) && (x[n])) {
m = x[n];
}
}
plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
mout = mxGetPr(plhs[0]);
*mout = m;
} else if (f) {
plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
mout = mxGetPr(plhs[0]);
*mout = m;
} else {
plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
}


}

END---cut here---


--
Dr Tristram J. Scott
Energy Consultant

tristram.scott@ntlworld.com (Tristram Scott) wrote in
message <UuOzj.25992$ab5.15398@newsfe1-win.ntli.net>...
> Tristram Scott <tristram.scott@ntlworld.com> wrote:
> > Steven Lord <slord@mathworks.com> wrote:
> >>
> >> "Tristram Scott" <tristram.scott@ntlworld.com> wrote
in message
> >> news:zjAzj.63245$os2.35954@newsfe3-win.ntli.net...
>
> [snip]
>
> >>
> >> That works for most cases. Being somewhat nitpicky,
though, I see that
> >> minnz(NaN) returns Inf and minnz([]) and minnz(0)
also return Inf. NaN is
> >> not equal to 0, so the code should return NaN in the
first case, and in
> >> the
> >> latter cases the result should be [], as the matrices
for those cases have
> >> no nonzero elements.
> >>
> >
> > Well Steve, it's nice to know that someone reads my
rubbish from time to
> > time. Certainly the code should check for those
trivial cases. Perhaps
> > I'll tackle that in the morning...
> >
>
> Well, the code has grown, but with little loss of
efficiency in the tests I
> have done.
>
> BEGIN---cut here---
> /*
> * minnz.c - Find minimum non-zero value in an array
> *
> * function m = minnz(x)
> *
> * Tristram Scott
> * 05/03/2008
> *
> * $Id: minnz.c,v 1.3 2008/03/06 09:02:55 tristram Exp $
> */
> #include <stdlib.h>
> #include <math.h>
> #include "mex.h"
> /* the gateway function */
> void mexFunction( int nlhs, mxArray *plhs[],
> int nrhs, const mxArray *prhs[] ) {
> const mxArray *tmp1;
> double m, *mout, *x;
> int f, n;
>
> if(nrhs!=1) {
> mexErrMsgTxt("One input required.");
> }
>
> tmp1 = prhs[0];
> n = mxGetM(tmp1) * mxGetN(tmp1);
> if (!n) {
> plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
> return;
> }
> x = mxGetPr(tmp1);
> f = 0; /* Flag */
> while (n--) {
> if (x[n]) {
> m = x[n];
> f = 1; /* Found a non-zero value */
> if (!isnan(x[n])) {
> break;
> }
> }
> }
> if (n>=0) {
> while (n--) {
> if ((x[n] < m) && (x[n])) {
> m = x[n];
> }
> }
> plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
> mout = mxGetPr(plhs[0]);
> *mout = m;
> } else if (f) {
> plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
> mout = mxGetPr(plhs[0]);
> *mout = m;
> } else {
> plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
> }
>
>
> }
>
> END---cut here---
>
>
> --
> Dr Tristram J. Scott
> Energy Consultant


Hi Tristram,

Thanks for preparing this new version.

yesterday i tried with m=infinity but it did not work so i
just replaced it by m=999999999999 (or so :-). like this i
was able to create the mex file. I tested it and i found
that minnz was almost 50% faster than using min(A(A>0))

today i tried with this new version of minnz and it does
reduce the time from 1528 to 884 seconds!!! :-)

Thanks a lot to you (and to everybody that helped me) !!!