Speed up recursive loop

I have the following loop which fills in a vector recursively. I need to speed it up as much as possible since it takes 60% of the computational time.

nobs    = 1e5;              % number of observations
p       = 2;                % number of lags
pseudoX = zeros(nobs,1);    % vector to fill in 
ro      = [0.1 0.5];        % autoregressive coefficients
noise   = rand(nobs,1)-p;   % white noise  

pseudoX(1:p) = randn(p,1);  % seed the vector
for ii = p+1:nobs
    pseudoX(ii) = ro * pseudoX(ii-p:ii-1) + noise(ii-p);
end

Maybe some mex help to get started (I work on win7 64bit and have MS Visual 2008).

EDIT Timings Win7 64bit 2012a-pre:

Solution to improve: Elapsed time is 0.233248 seconds.
Teja's filter      : Elapsed time is 0.003585 seconds.
Jan's decomposition: Elapsed time is 0.007797 seconds.

So far, Teja's solution is faster but Jan's is readable by a wider audience. I will wait the tests on previous releases to accept the answer and to see if I need flexibility on p.

[Jan]

Is p fixed?

[Oleg Komarov]

I might need some flexibility, but not at the moment. For now is fixed once I estimated the parameter.

Answer 1

[Jan]

Is this faster or slower:

ro1 = 0.1;
ro2 = 0.5;
for ii = p+1:nobs
  pseudoX(ii) = ro1 * pseudoX(ii-2) + ro2 * pseudoX(ii-1) + noise(ii-2);
end

[EDITED]: A C-mex - not optimized to allow for adjustments:

// Author: Jan Simon, Heidelberg, (C) 2012 [e-mail], BSD license
#  include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  double *ro, *noise, *seed,   // inputs
         ro1, ro2, k,          // parameters
         *Out;                 // output
  mwSize nobs, p, i;

    // Check number of inputs:
    if (nrhs != 5) {
       mexErrMsgIdAndTxt("JSimon:Func:BadNInput",
                         "5 inputs required.");
    }

    // Read inputs:
    nobs  = (mwSize) mxGetScalar(prhs[0]);
    p     = (mwSize) mxGetScalar(prhs[1]);
    ro    = mxGetPr(prhs[2]);
    noise = mxGetPr(prhs[3]);
    seed  = mxGetPr(prhs[4]);

    // Create output:
    plhs[0] = mxCreateDoubleMatrix(nobs, 1, mxREAL);
    Out     = mxGetPr(plhs[0]);

    // Copy seed:
    memcpy(Out, seed, p * sizeof(double));

    // The computations:
    ro1 = ro[0];
    ro2 = ro[1];
    k   = Out[p - 1];
    for (i = p; i < nobs; i++) {
       k      = ro1 * Out[i-2] + ro2 * k + noise[i-2];
       Out[i] = k;
    }

    return;
  }

This function is called as:

nobs  = 1e5;
p     = 2;
ro    = [0.1, 0.5];
noise = rand(nobs,1) - p;
seed  = randn(p, 1);

pseudoX = Func(nobs, p, ro, noise, seed)

On Win7, Matlab 2009a/64, MSVC 2008 this is 5 times faster than the FILTER method.

[Oleg Komarov]

Thanks Jan, 40x faster than my initial approach and one good example for me to start coding with mex.

Answer 2

[Teja Muppirala]

Do you know about the command "FILTER" from the Signal Processing Toolbox? This operation can be done very quickly.

nobs    = 1e5;              % number of observations
p       = 2;                % number of lags
pseudoX = zeros(nobs,1);    % vector to fill in 
ro      = [0.1 0.5];        % autoregressive coefficients
noise   = rand(nobs,1)-p;   % white noise  
pseudoX(1:p) = randn(p,1);  % seed the vector

tic
for ii = p+1:nobs
    pseudoX(ii) = ro * pseudoX(ii-p:ii-1) + noise(ii-p);
end
toc

tic
P2 = filter(1,[1 -ro(end:-1:1)],[pseudoX(1); pseudoX(2)-ro(2)*pseudoX(1); noise(1:end-2)]);
toc

max(abs(pseudoX-P2))

Pasting this into the editor and running it gives me:

Elapsed time is 0.618618 seconds.
Elapsed time is 0.003956 seconds.

ans =

       0

[Jan]

FILTER is part of the standard toolbox, so you do not require the SPT.

Unfortunately FILTER is not implemented efficiently (in Matlab <= 2011b). Hard-coded C will be much faster.

[Oleg Komarov]

Unfortunately I don't go further than a moving average with the filter function although I suspected it could have been used for such a problem.

Since I don't get completely what's going on, should concentrate on the matlab help to understand the algorithm or would you suggest some reference?

[Jan]

You can read the doc, of course. You find a Matlab implementation of FILTER here: http://www.mathworks.com/matlabcentral/answers/9900-use-filter-constants-to-hard-code-filter

A C-method which is faster than Matlab's filter in many cases: http://www.mathworks.com/matlabcentral/fileexchange/32261-filterm

[Oleg Komarov]

I tried to understand filter algorithm but I don't get how to generalize it.

[Oleg Komarov]

Thanks Teja, I really appreaciate your solution and although I think I understand how the algorithm works, I definitely have to spend some more time on the theory before feeling comfortable at using it.