From: Raphael on
"vortse a" <sonoffeanor-remove(a)yahoo.com> wrote in message <hq9lve$2gh$1(a)fred.mathworks.com>...
> Yi Cao has made a typo, he means to use the "diag" command. However there are two issues:
>
> 1. The results do not check. I don't know why.

so symmetrizing Neigh should make it work (Neigh = Neigh + Neigh').
From: Yi Cao on
"Raphael " <raphael.keller(a)gmail.com> wrote in message <hq9mt0$iba$1(a)fred.mathworks.com>...
> "vortse a" <sonoffeanor-remove(a)yahoo.com> wrote in message <hq9lve$2gh$1(a)fred.mathworks.com>...
> > Yi Cao has made a typo, he means to use the "diag" command. However there are two issues:
> >
> > 1. The results do not check. I don't know why.
>
> so symmetrizing Neigh should make it work (Neigh = Neigh + Neigh').

Yes, I assume Neigh is symmetric although the example I gave was not symmetric. For unsymmetric Neigh, it needs to be corrected as

AverageVelocity1 = diag(1./sum(Neigh,2))*Neigh*Velocity;

In terms of time, vectorization does not mean fast code any more due to the introducing of JIT (just-in-time) accelleration.

HTH
Yi
From: Bruno Luong on
"Raphael " <raphael.keller(a)gmail.com> wrote in message <hq99e6$736$1(a)fred.mathworks.com>...
> Hello,
>
> This is what I am trying to do : I have an n*3 matrix of cell velocities in space 'Velocity' and an n*n logical matrix of neighbourhood relationships 'Neigh' (0 if stranger, 1 if neighbour). n is about 5000 to 10000. Some speeds are NaN.
>
> To average speed in space, I did the following :
>
> AverageVelocity=NaN(n, 3);
> for i=1:n
> AverageVelocity(i,:)=nanmean(Velocity(Neigh(:,i),:));
> end

How many "1" Neigh has by column in average? If it's much smaller than n, you should take a look of SPARSE matrtix.

Bruno
From: Yi Cao on
"Yi Cao" <y.cao(a)cranfield.ac.uk> wrote in message <hq9tts$sp5$1(a)fred.mathworks.com>...
> "Raphael " <raphael.keller(a)gmail.com> wrote in message <hq9mt0$iba$1(a)fred.mathworks.com>...
> > "vortse a" <sonoffeanor-remove(a)yahoo.com> wrote in message <hq9lve$2gh$1(a)fred.mathworks.com>...
> > > Yi Cao has made a typo, he means to use the "diag" command. However there are two issues:
> > >
> > > 1. The results do not check. I don't know why.
> >
> > so symmetrizing Neigh should make it work (Neigh = Neigh + Neigh').
>
> Yes, I assume Neigh is symmetric although the example I gave was not symmetric. For unsymmetric Neigh, it needs to be corrected as
>
> AverageVelocity1 = diag(1./sum(Neigh,2))*Neigh*Velocity;
>
> In terms of time, vectorization does not mean fast code any more due to the introducing of JIT (just-in-time) accelleration.
>
> HTH
> Yi

The timing issue does not come from diag, but from calculation sequency. In the orginal code, I gave, it will firstly calculate diag(1./sum(Neigh,2))*Neigh, which is n^3 operations. Then, the results multiply Velocity, which is 3n^2 operations.

If we change the calculation sequence as follows:

AverageVelocity = diag(1./sum(Neigh,2))*(Neigh*Velocity);

It will first calculate (Neigh*Velocity), which is 3n^2 operations. Then the result will multiply by diag(1./sum(Neigh,2)), again is 3n^2.

Therefore, by changing calculation sequence, we can reduce operations from (n+3)n^2 to 6n^2. For large n, this will significantly reduce time. It should be compatible with or even faster than the loop based code.

HTH
Yi

From: Raphael on

> Therefore, by changing calculation sequence, we can reduce operations from (n+3)n^2 to 6n^2. For large n, this will significantly reduce time. It should be compatible with or even faster than the loop based code.
>
> HTH
> Yi
>

It's not faster, but no more than 10% slower.
Thank you very much for your answer.