Efficient Weighted Moving average

Question

Matlab User 2017-3-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329590-efficient-weighted-moving-average

编辑： Jan 2017-3-21

V.mat

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?

2 个评论
显示无隐藏无

José-Luis 2017-3-13

Could you provide example input and output?

Matlab User 2017-3-13

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 2017-3-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329590-efficient-weighted-moving-average#answer_259061

编辑：John BG 2017-3-16

在 MATLAB Online 中打开

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:jgb2012@sky.com jgb2012@sky.com>

9 个评论
显示 7更早的评论隐藏 7更早的评论

dpb 2017-3-19

编辑：dpb 2017-3-19

在 MATLAB Online 中打开

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>

John BG 2017-3-20

编辑：John BG 2017-3-20

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

what does addendum mean?

dpb 2017-3-20

编辑：dpb 2017-3-21

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 2017-3-13

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329590-efficient-weighted-moving-average#answer_258475

编辑：dpb 2017-3-17

在 MATLAB Online 中打开

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');

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2017-3-20

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

请先登录，再进行评论。

Answer 3

Jan 2017-3-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/329590-efficient-weighted-moving-average#answer_259567

编辑：Jan 2017-3-20

在 MATLAB Online 中打开

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.

2 个评论
显示无隐藏无

dpb 2017-3-20

编辑：dpb 2017-3-21

在 MATLAB Online 中打开

"...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 2017-3-21

编辑：Jan 2017-3-21

在 MATLAB Online 中打开

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.

请先登录，再进行评论。