converting for loop to matrix operation

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

[Sara]

You'll need to post the function daily_harmonic_regression to get help

[ahmed shaaban]

I attached the code

if true
function [anom,cycle]=daily_harmonic_regression(data)
t=(1:length(data))';
sint=sin(2*pi*t/365.25);
cost=cos(2*pi*t/365.25);
sin2t=sin(2*pi*2*t/365.25);
cos2t=cos(2*pi*2*t/365.25);
sin3t=sin(2*pi*3*t/365.25);
cos3t=cos(2*pi*3*t/365.25);
sin4t=sin(2*pi*4*t/365.25);
cos4t=cos(2*pi*4*t/365.25);
X=[ones(length(data),1) sint cost sin2t cos2t sin3t cos3t sin4t cos4t];
C=inv(X'*X)*X'*data;
cycle=X*C;
anom=data-cycle;
end

[dpb]

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]

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]

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);

[ahmed shaaban]

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]

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!!!