converting for loop to matrix operation

Question

ahmed shaaban 2014-5-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/131656-converting-for-loop-to-matrix-operation

评论： Andrei Bobrov 2014-6-2

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 个评论
显示 2更早的评论隐藏 2更早的评论

dpb 2014-5-29

在 MATLAB Online 中打开

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 2014-6-2

Thanks a lot. You mean that after calculating C, I need to write solver for 1 dimension, 2 dimensions , and so on .

请先登录，再进行评论。

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

Answer 1

Andrei Bobrov 2014-5-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/131656-converting-for-loop-to-matrix-operation#answer_138827

编辑：Andrei Bobrov 2014-6-2

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

ahmed shaaban 2014-6-2

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.