bugGNU Octave - Bugs: bug #65245, movfun/movmean: scope for speed up...

 
 

bug #65245: movfun/movmean: scope for speed up with conv

Submitter:  A.R. Burgers <arb>
Submitted:  Sat 03 Feb 2024 12:18:16 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Feature Request
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Fri 15 Mar 2024 08:02:01 PM UTC, comment #5: 

Can this then be assigned to somebody?
We could request a patch from the author of the original post.

It's been a long time since I modify Octaves core code, so I wouldn't feel comfortable doing this myself. Also, I still have limited resources and would react very slowly.
If nobody else can do it, I'll take the shot.

I am changing the Item Group to "Feature Request".

Juan Pablo Carbajal <juanpi>
Group Member
Fri 15 Mar 2024 07:02:43 PM UTC, comment #4: 

(we don't currently support the "property=value" syntax)

but since the current compatible form of movfun includes

Y = movfun (..., "PROPERTY", VALUE)

it would be low risk of a compatibility collision by adding such an option as a octave-unique property.

Nicholas Jankowski <nrjank>
Group Member
Fri 15 Mar 2024 06:34:33 PM UTC, comment #3: 

One alternative, given that the mfile implementation is quite straight forward is to provide a "algorithm" or "method" to choose how to do this calculations.
The default value for those arguments would be "movfun", but the movxxx functions that yield to the conv method, could accept "conv" as a value for that argument.

So movmean(a, wlen) is as of now, but movmena(a, wlen, method="conv") would use the proposed method.

JPi

Anonymous
Fri 09 Feb 2024 02:31:07 PM UTC, comment #2: 

No, not so easy to integrate this in the current framework, just wanted to have a record of the idea for the time being.

movstd and movvar also yield to this approach


np = 100;
nw = 13;
window = ones(1, nw);
a1 = rand(1, np);
a2 = a1 .^2;
a0 = ones(1, np);
count  = conv(a0, window, 'same');
count1 = count - 1;

fc = @(x) conv(x, window, 'same');

conv_sum  = fc(a1);
conv_mean = fc(a1) ./ count;
conv_var0 = (fc(a2) - conv_mean .* (2 * fc(a1)  - nw * conv_mean)) ./ count1;
conv_var1 = (fc(a2) - conv_mean .* (2 * fc(a1)  - nw * conv_mean)) ./ count;
conv_std0 = sqrt(conv_var0);
conv_std1 = sqrt(conv_var1);

mov_sum  = movsum(a1, nw);
mov_mean = movmean(a1, nw);
mov_var0 = movvar(a1, nw);
mov_var1 = movvar(a1, nw, 1);
mov_std0 = movstd(a1, nw);
mov_std1 = movstd(a1, nw, 1);

assert(all(abs(mov_sum - conv_sum)) < 1.e-10);
assert(all(abs(mov_mean - conv_mean)) < 1.e-10);
assert(all(abs(mov_var0 - conv_var0)) < 1.e-10);
assert(all(abs(mov_var1 - conv_var1)) < 1.e-10);
assert(all(abs(mov_std0 - conv_std0)) < 1.e-10);
assert(all(abs(mov_std1 - conv_std1)) < 1.e-10);


A.R. Burgers <arb>
Sat 03 Feb 2024 06:31:33 PM UTC, comment #1: 

Yes, this is true.  It might be a good addition, but there is also a bigger strategic issue that is worth thinking about.  The current architecture is that the movXXX functions are simple in that they do just a little bit of input validation and then pass on the function to movfun().  This structure makes maintenance easy and any improvements to movfun() are automatically reflected in the derivative movXXX functions.  movmean() is a special case in that there is a convenient kernel that can be used with conv() (which is written in C++) to accelerate that function.  But other functions, such as movmad(), do not seem to have any easy mapping on to conv().  Perhaps it would be better to re-write movfun() in C++ where loops would not be a problem.

Rik <rik5>
Group administrator
Sat 03 Feb 2024 12:18:16 PM UTC, original submission:  

see also bug #55389

with a a vector, window size wlen odd, conv can give much faster results because of the linearity of the operator in the movsum/movmean cases.

The following two snippets show how conv can calculate the same result for movsum and movmean


result = movsum(a, wlen);
result = conv(a, ones(1, wlen), 'same');



result = movmean(a, wlen)
count = conv(ones(size(a)), ones(1, wlen), 'same');
result = conv(a, ones(1, wlen), 'same') ./ count;


Speed-ups can be a factor ~30


movsum time 0.321551
movsum with conv time 0.0108559
movmean time 0.520523
movmean with conv time 0.0155551


Timings obtained with this script


np = 10000; ns = 13;
nrep = 100;
a = rand(1, np);
a1 = ones(1, np);

tic;
for i = 1 : nrep
  as = movsum(a, ns);
end
fprintf('movsum time %g\n', toc);

tic;
for i = 1 : nrep
  as_c = conv(a, ones(1, ns), 'same');
end
fprintf('movsum with conv time %g\n', toc);

assert(all(abs(as - as_c)) < 1.e-10);

tic;
for i = 1 : nrep
  as = movmean(a, ns);
end
fprintf('movmean time %g\n', toc);

tic;
for i = 1 : nrep
  as_c = conv(a, ones(1, ns), 'same') ./ conv(a1, ones(1, ns), 'same');
end
fprintf('movmean with conv time %g\n', toc);


A.R. Burgers <arb>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by arb (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-03-15 juanpi Item GroupPerformance Feature Request
    2024-02-03 rik5 StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-4b48.
    Corresponding source code