Speed up recursive loop
9 views (last 30 days)
Show older comments
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.
Accepted Answer
Jan
on 8 Feb 2012
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 j@n-simon.de, 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.
More Answers (1)
Teja Muppirala
on 8 Feb 2012
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
5 Comments
See Also
Categories
Find more on Whos in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!