Vectorized array operation which depends on previous array value

Hi,

There is any way of making this operation vectorized?

a = ones(1000,1);
b = rand(1000,1);
c = rand(1000,1);
for i=2:1000
    a(i) = b(i) + a(i-1) .* c(i);
end

I really appreciate your help, Many thanks, Dylan

[Matt J]

I think you would have to use either Simulink or a MEX routine.

Answer 1

[Jan]

Do you have a C compiler installed? Then a C-mex would be the best option.

// File: YourFcn.c
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
  size_t n, i;
  double a1, *b, *c, *r;

    n  = mxGetNumberOfElements(prhs[1]);
    a1 = mxGetScalar(prhs[0]);
    b  = mxGetPr(prhs[1]);
    c  = mxGetPr(prhs[2]);

    plhs[0] = mxCreateDoubleMatrix(n, 1, mxREAL);
    r       = mxGetPr(plhs[0]);

    r[0]    = a1;
    for (i = 1; i < n; i++) {
      r[i] = b[i] + r[i-1] * c[i];
    }
  }

UNTESTED! Running this function might crash your computer. Call it like:

a = YourFcn(a1, b, c);

[EDITED] Now the first element of a is provided as input.

For productive usage, add some checks of the inputs arguments: if the numbers of elements of b and c are equal, is the inputs are doubles and do not have an imaginary part, if they are sparse arrays, ... These checks are very cheap in a Mex function, but forgetting them might cost you hours for debugging.

[Dylan Marques]

Thanks Jan, it seems the better option to me

[Walter Roberson]

a(1) is needed for input. The value can make quite a bit of difference to the final result.

[Jan]

@Walter: In the original code a is set to ones completely, but actually only a(1) matters. I have added the input of a(1) now instead of setting it to 1.0 .

Answer 2

[Torsten]

a(2:end)=b(2:end)+a(1:end-1).*c(2:end);

Best wishes

Torsten.

[Walter Roberson]

This would copy out the original a values and use them in the computation. The original code uses the a values that were just set in the previous iteration.

[Dylan Marques]

Yeah, I need sort of a cumsum for a personalized function.