Efficient Weighted Moving average

I have a matrix with x, y, and an uncertainty value (25x3 double). I want to compute a moving average of these points weighted by the third value, (V(:,3)), with a time window of say 5 points. What is the most efficient way to achieve this?

[José-Luis]

Could you provide example input and output?

[Matlab User]

I have attached the input as V=x,y,uncertainty. The uncertainty is what I wish to weight the average by - the lower the value the higher the weight.

Answer 1

[John BG]

Hi Matlab User

1.

simulating data

clear all,clc

V=zeros(25,3);x=V(:,1);y=V(:,2);
w=randi([1 9],25,1);
w=w/sum(w);V(:,3)=w

V(:,1)=randi([-10 10],25,1)   % x
V(:,2)=randi([-10 10],25,1)   % y
V(:,3)=w;

.

the data

format short

V =
   -3.0000    3.0000    0.0857
   -5.0000    8.0000    0.0190
  -10.0000    3.0000    0.0476
   -2.0000   -1.0000    0.0476
    4.0000    1.0000    0.0571
         0    3.0000    0.0190
    1.0000   -7.0000    0.0190
   10.0000   -1.0000    0.0857
   10.0000   -1.0000    0.0762
   -6.0000   -3.0000    0.0762
   -5.0000   -9.0000    0.0381
         0    4.0000    0.0095
   -7.0000    1.0000    0.0190
   -8.0000  -10.0000    0.0286
    3.0000    8.0000    0.0381
   -4.0000   -4.0000    0.0095
    2.0000    6.0000    0.0190
   -9.0000         0    0.0095
    5.0000    6.0000    0.0286
   -4.0000   -2.0000    0.0286
  -10.0000    6.0000    0.0667
   -9.0000   -2.0000    0.0381
   -9.0000    1.0000    0.0095
    7.0000    4.0000    0.0476
    9.0000   -6.0000    0.0762

2.

averaging with the shifting window and applying the weight V(:,:,3)

N=25
R=zeros(N,1)
for k=3:1:N-2
    s=[k-2 k-1 k k+1 k+2]
    R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3))
end

R =
         0
         0
   -0.0133
   -0.0076
   -0.0248
    0.0857
    0.1686
    0.0714
    0.0124
    0.0276
   -0.0610
   -0.1810
   -0.0705
   -0.0248
   -0.0133
   -0.0105
    0.0724
    0.0133
   -0.0057
   -0.0629
   -0.0619
   -0.0410
   -0.0010
         0
         0

This does not apply the shifting window on the first 2 and last 2 sames, there aren't 5 samples to apply the window on head and tail anyway, aren't they?

if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?

To any other reader, please if you find this answer

please click on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[Matlab User]

That's great - thankyou for your time!

[dpb]

>> N=5;wt=ones(1,N)/N;                           % compute MA weights
>> C=conv(mean(V(:,1:2),2).*V(:,3),wt,'valid');  % N-wt MA of Avg.*WGT
>> Z=zeros(fix(N/2),1); [[Z;C;Z] R]              % compare looping sol'n
ans =
       0         0
       0         0
 -0.0134   -0.0133
 -0.0077   -0.0076
 -0.0248   -0.0248
  0.0857    0.0857
 ...        ...
 -0.0619   -0.0619
 -0.0410   -0.0410
 -0.0010   -0.0010
       0         0
       0         0
>>

Use the vectors, luke! (OP did ask for 'efficient', after all... ) :)

[John BG]

Hey dpd

I am trying to tic toc your answer against mine, with same example signal I generate above, but there's this horzcat error message after you answer

would you please be so kind to clarify?

    ..
    tic
    R=zeros(N,1);
    tic
    for k=3:1:N-2
        s=[k-2 k-1 k k+1 k+2];
        R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3));
    end
    toc
    tic
    wt=ones(1,N)/N;
    C=conv(mean(V(:,1:2),2).*V(:,3),wt,'valid');  % N-wt MA of Avg.*WGT
    Z=zeros(fix(N/2),1); [[Z;C;Z] R] 
    toc
Error using horzcat
Dimensions of matrices being concatenated are not consistent.

.

John BG

<mailto:[e-mail] [e-mail]>

[John BG]

DPD: also, variable R not defined

John BG

<mailto:[e-mail] [e-mail]>

[Jan]

@John BG: The R is the result of your code. Omit the "Z=zeros(fix(N/2),1); [[Z;C;Z] R]" line, because it is for comparing the results only. See dpb's answer.

[dpb]

Jan 'splained the R issue; if your result is the full 25 elements there should be no issue. But, as Jan notes, that's only there to compare the results side-by-side visually, no bearing on anything else.

>> tic;N=25;
R=zeros(N,1);
for k=3:1:N-2
  s=[k-2 k-1 k k+1 k+2];
  R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3));
end;toc
Elapsed time is 0.019355 seconds.
>> tic;C=conv(mean(V(:,1:2),2).*V(:,3),wt,'valid');toc
Elapsed time is 0.000394 seconds.
>> 
>> 0.019355/0.000394
ans =
 49.1244
>>

on a 15-yo 32-bit dinosaur...

JIT will undoubtedly help some if put in m-files, but the builtins will still win hands down I'm sure...

ADDENDUM

Oh, that's not quite fair, I just retrieved the compute line w/o the overhead...let's see if make more ephen-stephen footing--

>> tic;N=5;C=conv(mean(V(:,1:2),2).*V(:,3),ones(1,N)/N,'valid');toc
Elapsed time is 0.001087 seconds.
>>

Ah, that's more real--

>> 0.019355/0.001087
ans =
 17.8059
>>

Still sizable advantage altho ONES is a killer it seems...hmmm, another test or two and at least at command line [1 1 1 1 1]/5 doesn't help much, either. I knew '[]' are considered harmful to speed from the editor optimization hints; wasn't aware were such a killer.

I didn't put into m-files, however; will leave that as "exercise for Student" :)

ADDENDUM 2:

It might be interesting if you want to delve into where issues are to time only the loop just to see if by any chance't the overhead is the issue. I don't think it will be, but "ya' never know unless ya' asks"...

Also, might try some minor efficiencies in the loop and see if can help it out any--

N1=fix(N/2);     % where to start loop
s=[N1-1:N1+1];   % initial indexing vector
for k=N1+1:N-2
  R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3));
  s=s+1;         % just increment s instead of rebuilding every loop...
end

You can also try rewriting the calculation as well...

R(k)= mean(mean(V(s,1:2),2)).*V(s,3));

Probably some other ways to help a tad as well, but the key item will still be found to avoid the loop entirely I suspect, using instead the builtin routines that are compiled wherever can...but, it's a learning experience in becoming adept using Matlab as it is intended everybody has to work thru initially. Of course, it helps to have good background in matrix algebra, numerical analysis, signal processing, etc., etc., etc., .... :)

A quote at least attributed to if not verified as actually his from Mark Twain may be apropos--

"Good decisions come from experience. Experience comes from making bad decisions." <vbg>

[dpb]

PS: @John BG, my initials are dp*B*, not dpd...just for the record ;)

[John BG]

ok DPB, may I use capitals, or do you want them small, in small characters only?

what does addendum mean?

[dpb]

However is fine... :) I've noticed quite a few have typed 'dpd' rather than the trailing 'b'; I'm not sure why. I kinda' like the geometry of lower case string however as being mirror image/rotations of a single basic form...

"Addendum. An addendum, in general, is an addition required to be made to a document by its author subsequent to its printing or publication. It comes from the Latin verbal phrase addendum est, being the gerundive form of the verb addo, addere, addidi, additum, "to give to, add to", meaning "(that which) must be added"."

Answer 2

[dpb]

N=5;
yavN=conv(v(:,2).*v(:,3),ones(1,N)/N,'valid');

discarding the zero-padded end values. See

doc conv  % for other options, details...

ADDENDUM

If not clear how to use mean in the above, simply compute it "on the fly" as well instead of using just a single column from V...note the use of the optional DIM second argument to the MEAN function...

C=conv(mean(V(:,1:2),2).*V(:,3),ones(1,N)/N,'valid');

[Jan]

+1. This is much faster than the methods using loops.

Answer 3

[Jan]

Some further speed measurements:

function speedtest
N = 10000;
tic; for k = 1:N; R = test1(V); end, toc
tic; for k = 1:N; R = test2(V); end, toc
tic; for k = 1:N; R = test3(V); end, toc
tic; for k = 1:N; R = test4(V); end, toc

function R = test1(V)  % John BG 
N=size(V, 1);
R=zeros(N,1);
for k=3:1:N-2
   s=[k-2 k-1 k k+1 k+2];
   R(k)= mean(.5*(V(s,2)+V(s,1)).*V(s,3));
end

function R = test2(V)  % Jan
% 1. repeated calcultations moved outside the loop
% 2. Indexing by a:b is much faster than using a pre-calculated vector
N = size(V, 1);
R = zeros(N,1);
X = 0.5 * sum(V(:, 1:2), 2) .* V(:, 3);
for k = 3:1:N-2
   R(k) = sum(X(k-2:k+2)) * 0.2;
end

function R = test3(V)   % dpb
N = 5;
R = conv(mean(V(:,1:2), 2) .* V(:,3), ones(1,N) / N,'valid');

function R = test4(V)   % dpb without MEAN()
N = 5;
R = conv(0.5 * (V(:,1) + V(:, 2)) .* V(:,3), ones(1,N) / N, 'valid');

Results (R2016b/64, Core2Duo):

 Elapsed time is 3.138929 seconds.   John BG
 Elapsed time is 0.444694 seconds.   Improved loop
 Elapsed time is 0.193468 seconds.   dpb
 Elapsed time is 0.119400 seconds.   dpb without MEAN()

What happens, if the input is not tiny:

V = rand(25000, 3);
N = 10;

 Elapsed time is 3.667920 seconds.
 Elapsed time is 0.418300 seconds.
 Elapsed time is 0.018349 seconds.
 Elapsed time is 0.014222 seconds.

Conclusion: Loops profit from smart indexing and avoiding repeated calculations, but conv rules. mean has a remarkable overhead.

[dpb]

"...mean has a remarkable overhead."_

Interesting, wonder why so much. Well, let's see...

C:\ML_R2014b\toolbox\matlab\datafun\mean.m
>> edit mean

Huh. It's not a builtin that I'd just presumed would be given the ubiquitousness of the operation. Has lots of input checking as well before finally getting around to doing any work.

Have to remember that when speed is of essence (altho since I've given up the consulting gig for the farm, that's not likely to ever be so for me again, but for others...)

ADDENDUM

Wonder how much penalty there would be in

sum(V(:,1:2),2)/2

in comparison as replacement syntax for explicit sum (I hate having to write V twice't when subscripting can be used but that's just me, granted)...

[Jan]

x = rand(1000, 2);
tic; for k = 1:100000; m = x(:, 1) + x(:, 2); end; toc
tic; for k = 1:100000; m = sum(x, 2);         end; toc

 Elapsed time is 0.120969 seconds.
 Elapsed time is 0.156458 seconds.

The difference can be measured. SUM is parallelized for more than 88999 elements. Then the number of cores matter also.

x = rand(100000, 2);
tic; for k = 1:1000; m = x(:, 1) + x(:, 2); end; toc
tic; for k = 1:1000; m = sum(x, 2); end; toc

 Elapsed time is 0.103634 seconds.
 Elapsed time is 0.148599 seconds.

But rmember, we are talking about 0.00004 seconds. "sum(x,2)" can save minute of hourse during the debugging. Therefore I prefer sum and apply the dirty tricks like v(:,1)+v(:,2) only for bottlenecks and with an exhaustive documentation and a unit-test function, which specifically cares about these lines.