converting for loop to matrix operation
8 views (last 30 days)
Show older comments
ahmed shaaban
on 29 May 2014
Commented: Andrei Bobrov
on 2 Jun 2014
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
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.
Accepted Answer
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
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!!!
More Answers (0)
See Also
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!