converting for loop to matrix operation

10 views (last 30 days)
Hello All I wondering if I can do the following calculation without using for loop I am using function called daily_harmonic_regression which takes one vector data, that the reason that I use for loop. For big matrix, such calculation takes a lot of time. msl is (length(long),length(lat), length(time))
anom_mslp=msl;
for lo=1:length(long);
for la=1:length(lat)
[anom,cycle_mslp]=daily_harmonic_regression(squeeze(msl(lo,la,:)));
anom_mslp(lo,la,:)=anom;
end
end
thanks in advance
  4 Comments
dpb
dpb on 29 May 2014
function [anom,cycle]=daily_harmonic_regression(data)
t=2*pi*(1:length(data))/365.25.';
sint=sin(t);
cost=cos(t);
sin2t=sin(2*t);
etc., ...
X=[ones(length(data),1) sint cost sin2t cos2t sin3t cos3t sin4t cos4t];
C=inv(X'*X)*X';
Up to here for a given array everything is redundant for every column so factor this out into another routine that returns C. Then you can write a second solver routine the applies it to each column of data.
BTW, the inverse of the normal equations is a pretty numerically marginal solution technique; may want to look at Matlab '\' operator.
ahmed shaaban
ahmed shaaban on 2 Jun 2014
Thanks a lot. You mean that after calculating C, I need to write solver for 1 dimension, 2 dimensions , and so on .

Sign in to comment.

Accepted Answer

Andrei Bobrov
Andrei Bobrov on 30 May 2014
Edited: Andrei Bobrov on 2 Jun 2014
function [anom,cycle]=dhr(data) % daily_harmonic_regression
d = data(:);
n = numel(d);
a = 2*pi*(1:n)'*(1:4)/365.25;
ii = reshape(1:8,[],2)';
x0 = [sin(a), cos(a)];
x = [ones(n,1), x0(:,ii(:))];
cycle = reshape(x*((x'*x)\x'*d),size(data));
%{
(x'*x)\x' or inv(x'*x)*x' -> matrix is close
to singular or badly scaled for any data - always!!!
%}
anom=data-cycle;
end
use dhr
[z,y] = cellfun(@dhr,num2cell(msl,3),'un',0);
anom_mslp = cell2mat(z);
ADD
function [anom,cycle]=dhr0(data) % daily_harmonic_regression
[m,n,k] = size(data);
a = 2*pi*(1:k)'*(1:4)/365.25;
i0 = reshape(1:8,[],2)';
x0 = [sin(a), cos(a)];
x = [ones(k,1), x0(:,i0(:))];
X0 = (x'*x)\x';
cycle = zeros(m,n,k);
anom = zeros(m,n,k);
for ii = 1:m
for jj = 1:n
d = data(ii,jj,:);
cycle(ii,jj,:) = x*(X0*d(:));
anom(ii,jj,:) = d - cycle(ii,jj,:);
end
end
%{
(x'*x)\x' or inv(x'*x)*x' -> matrix is close
to singular or badly scaled for any data - always!!!
%}
end
use dhr0
[anom,cycle]=dhr0(msl);
  2 Comments
ahmed shaaban
ahmed shaaban on 2 Jun 2014
Thanks a lot, it look like complicated to write the code in matrix form. I have two issues, first, the execution time for the code using iterations took 38.455499 seconds, and for your code it took 43.352565 seconds. second, there is some difference between output from both codes. I have attached the output for both. blue one is for your code. Thanks a lot.
Andrei Bobrov
Andrei Bobrov on 2 Jun 2014
First,please use function dhr0.m, see code in my answer after word ADD.
Second: (x'*x)\x' or inv(x'*x)*x' -> matrix is close to singular or badly scaled for any data - always!!!

Sign in to comment.

More Answers (0)

Categories

Find more on Creating and Concatenating Matrices 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!