Speed up recursive loop

2 次查看（过去 30 天）
Oleg Komarov
Oleg Komarov 2012-2-8
编辑： Jan 2023-1-29
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.
  2 个评论
Jan
Jan 2012-2-8
Is p fixed?
Oleg Komarov
Oleg Komarov 2012-2-8
I might need some flexibility, but not at the moment. For now is fixed once I estimated the parameter.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2012-2-8
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.
  1 个评论
Oleg Komarov
Oleg Komarov 2012-2-9
Thanks Jan, 40x faster than my initial approach and one good example for me to start coding with mex.

请先登录，再进行评论。

更多回答（1 个）

Teja Muppirala
Teja Muppirala 2012-2-8
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 个评论
显示 3更早的评论
Oleg Komarov
Oleg Komarov 2012-2-8
I tried to understand filter algorithm but I don't get how to generalize it.
Oleg Komarov
Oleg Komarov 2012-2-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by